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ABSTRACT 


Thia report describes the progress under NASA Cooperative 
Agreement NCCS-8 during the period from May I, 1978 to April 
30, 1979. The seismic wave velocity structure in the crust 
and upper mantle region beneath the Tibetan plateau was studied 
in detail. Also, a preliminary study of the uppermost mantle 
P wave velicity beneath Iran and Turkey was carried out, and 
the results are compared with those for the Tibetan plateau. 
These two studies compose the bulk of our efforts on the ob- 
servational aspects of continental collision zones in addition 
to satellite derived data. On the theoretical aspects the 
thermal evolution of converging plate boundaries was explored 
using a finite difference scheme. A separate study examined 
the question of whether continental lithoshpere can be suh 
ducted into the asthenosphere by a large amount. 
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INTRODUCTION 

A geophysical study of the structure and processes of 
the continental convergence zones with special references 
to the Alpine-Hixnalayan Belt is carried out under NASA Coop- 
erative Agreement NCC5-8. This report presents the accom- 
plishments made during the period of May 1, 1978 to April 
30, 1979. 

The bulk of this report consists of three published papers 
and one paper in press. The first two papers concentrate on 
obtaining reliable seismic data to complement satellite 
derived data. In the first paper, a combination of surface 
wave dispersion, (a leaky mode) propagation, refraction 
profiles using earthquakes, and teleseismic P and S wave 
travel time residuals place constraints on the seismic wave 
velocity structures beneath the Tibetan plateau. The results 
confirm the idea of a very thick crust ('^7 0 km) and reveal a 
shield- like uppermost mantle (with seismic velocities of 
8.1km/sfor P wave and 4.8 km/s for S wave) overlying a low 
velocity mantle. The uppermost mantle P wave velocities 
beneath Turkey and Iran are investigated in the second paper. 
These velocities show large regional variations (7.7 km/s 
for Turkey and 8.0 km/s for Iran) - an indication of possible 
differences in tectonic processes. These velocity values were 
further used to estimate the relative temperature near the Moho 
between these regions and stable shields. 

We also finished two theoretically oriented studies. The 
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evolution of the thermal structure beneath a subduction zone 
is modeled using numerical methods. This is a starting point 
of modeling collisional tectonics because, at least along the 
Alpine-Himalayan belt, subduction of oceanic lithospheres 
always proceded continental collision. Thus, the thermal 
regime of a subduction zone serves as the initial condition 
of the complicated thermal mechanical processes of a collision 
zone. The problem of whether the continental ' Ithosphere, 
with its buoyant crust, can be subducted by a large amount 
into the asthenosphere is investigated in the last study 
presented in this report. The uncertainties are large so that 
one can not prove this is definitely not the case. However, 
the most favorable condition for subducting continental 
lithosphere would involve detachunent of the continental 
crust from the rest of the lithosphere. 

Reprints of the published papers and preprint of the 
paper in press follow this section. In the end of this report 
we also include abstracts of two studies presented at pro- 
fessional meetings during this period. 


CONSTRAINTS ON THE SEISMIC WAVE 
VELOCITY STRUCTURE BENEATH THE TIBETAN PLATEAU 
AND THEIR TECTONIC IMPLICATIONS 
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Massachusetti: Institute of Technology 


Cambridge, MA 
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Abstract 

W« combina obsarvatlons of group and phase velocity 
dispersion of Rayleigh waves, of the wave form of a long 
period phase, of and velocities from unreversed 
refraction profiles using earthquakes, and of teleseismic 
S-P travel time resiv.ualj to place bounds on the seismic 
wave velocity structure of the crust and upper mantle under 
Tibet. From surface wave measurements alone, the Tibetan 
crustal thickness can be from 55 km to 85 km, with 
corresponding uppermost mantle shear wave velocities of 
about 4.4 km/s to 4.9 km/s, respectively. The P^ and 
velocities were determined to be 8.12;f0.06 km/s and 4.8;^0.1 
km/s respectively using travel time data at Lhasa from 
earthquakes in and on the margins of Tibet. With these 
estimates of the uppermost mantle velocity and with the 
surface wave dispersion, the crustal thickness is most 
likely to be between 65-80 km with an average shear wave 
velocity in the crust less than 3.5 km/s. A synthesis of 
one wave form does not provide an additional constraint 
on the velocity structure but is compatible with the range 
of models given above. Measurements oZ both teleseismic S 
wave and P wave arrival times for nine earthquakes within 
Tibet show unusually large intervals between P and S 
compared with the Jef f reys-Bullen Tables. Yet both S wave 
and P wave arrival times from eight earthquakes in the 
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nearby Himalaya agree quite well with these Tables. Thus 
the P,) and velocities apparently do not reflect high 
velocities in the mantle to a great depth beneath Tibet. 

Prom the dependence of the seismic velocities of olivine on 
pressure and temperature and frmt the similarity of the 
measured P^ and velocities beneath Tibet and beneath 
shields and platformr, the temperature at the Moho beneath 
Tibet is compatible with being 250-300** higher than beneath 
shields and platforms, i.e. 750*C if the temperature beneath 
the platforms is close to 500*C. Such a temperature could 
reach or exceed the solidus of the lower crust. Simple one- 
dimensional heat conduction calculations suggest that the 
tectonic and volcanic activity could be explained by the 
recovery of the geotherm maintained by a mantle heat flux 
of about 0.9 HFU at the base of the crust. If the 
distribution of radioactive heat production elements were 
not concentrated at the top of the crust, radioactive heating 
could also contribute significantly to the recovery of the 
geotherm and thur lower the required mantle heat flow. 

Thus the idea of a thickened crust in response to horizontal 
shortening is compatible both with these data and with these 
calculations . 
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INTRODUCTION 

With an average elevation of about S km above sea 
level over an area of 7 x 10'^ km‘, the Tibetan plateau is 
certainly one of the most conspicuous topographic '’matures 
on earth. To the south and southwest, it is bounded by the 
Himalayan convergent zone, where the continental collision 
between India and Eurasia plates manifests itself 
predominantly by low angle thrust faulting [e.g. Fitch, 1970; 
Molnar et al., 1973, 1977]. The Altyn Tagh and the Kun Lun 
left-lateral strike-slip faults form its northern boundary 
[e.g. Molnar and Tapponnier, 1975; Tapponnier and Molnar, 
1977], In contrast, recent normal faulting [e.g. Moln<.r 
and Tapponnier, 1975, 1978; Ni and York, 1978] and volcanism 
are widespread over the plateau [Burke et al., 1974). With 
its complex variations of active tectonic styles and its 
massive volume, the evolution of the Tibetan p'.ateau has 
been recognized as a key to the understanding of continent- 
continent convergent processes [e.g. Dewey and Burke, 1973; 
Molnar and Tapponnier, 1975, 1978; Toksflz and Bird, 1977]. 

Primarilv because of its inaccessibility, geological 
data available about Tibet are very limited. Upper Cretaceous 
limestones seem to cover much of southern Tibet [Hennig, 

1915; Norin, 1946]. This observation places a lower bound 
on the time of uplift, but we are not aware of evidence that 
constrains the timing more tightly. Volcanism, which is 
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app«r«ntly quit* young Is wldsspi'ssd on th« plstsau [s.g. 

Burko, «t si., 1974} Kidd, 1975]. Rscsntly Dsng [1978] 
rsportsd Qustsrnary calc-alkallna to alkalina volcanics in 
northarn Tlbat at sltas close to those of the earlier findings 
of Norin [1946]. Similar rock types have also been reported 
along the southern part of Tibet (Chang and Cheng, 1973; 

Hennig, 1915; Kidd, 1975]. 

Given the unusual elevation and the limited geologic 
data, the crustal and upper mantle structures should place 
important constraints on the geologic evolution of the plateau. 
Gravity data are too limited to provide a tight constraint. 
Amboldt (19481 reported one gravity measurement within the 
plateau that is consistent with isoscatic equilibrium. 

Chang and Cheng (1973] also inferred isostatic equilibrium 
from more recent, but unpublished gravity data. As pointed 
out by Bird and Toksdz (1977] the crustal thickness of Tibet 
is likely to be from 55 to 70 km assuming isostatic 
c^pensation with mantle properties varying from a mid-ocean 
ridge-like conditions to stable shield-like environments. 

Given the uncertainties in the density of the crustal rocks, 
an even greater range is possible. 

Almost all the previous seismic studies of Tibet are 
based on group velocity dispersion curves of surface waves 
(Bird, 1976; Bird and Toksdz, 1977; Chun and Yoshii, 1977; 
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Patton/ 1978; Tung and Tang, 1974; Tsenq and Sung/ 1963]. 

In ganaral/ thaa* invastigatora concludad that tha cruat 
ia about 70 km thick with a ralativaly low uppar mantla 
ahaar wava valoclty (naar 4.5 km/a). Tha ma.itla valocitiaa/ 
howavar/ cannot ba wall conatrainad with aurfaca wava data 
alona [Dar at al./ 1970). Tha poaaibla trada-off batwaan 
uppar mantla valocitiaa and cruatal thicknaaa waa not 
invaatigatad in moat of thaaa atudiaa/ and littla attantion 
haa baan givan to tha uncartaintiaa and bounda of allowabla 
velocity atructuraa. Wa do not conaidar tha raaulta in 
the publiahed atudiaa to ba concluaive avidanca for cruatal 
thickneaaaa of 70 kir>, and conaaquantly wa have carriad out 
a atudy of not only aurface wava diaparaion but alao of 
othar aaiamic phaaea. 

Conaidar ing that the cruat ia thick now, but probably 
waa not thick before tha collision between India and Eurasia 
in the early Tertiary, there have been two extreme mechanisma 
proposed for the thickening of the crust and the uplift of 
the plateau. Dewey and Burke [1973], Toksdz and Bird [1977] 
and probably others suggested that Tibet formed by relatively 
uniform crustal shortening and thickening (like an accordion) 
in response to horizontal compression. Others [e.g., Argand, 
1924; Powell and Conaghan, 1973, 19751 inferred that low 
angle underthrusting of one crustal block (India) beneath 
another (Tibet) is the dominant mechanism for creating the 
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high plAtMU. Both modols require e decrease » if only 
temporary, xn the crustal temperature gradient due to the 
crustal thickening. Therefore a more quantitative treatment 
of either situation is required to explain how volcanism 
could occur on the plateau. 

The existence oi the volcanoes and the inference of 
relatively low seismic wave velocities led logically to the 
suggestion that the upper mantle beneath Tibet is relatively 
hot and light. Part of the high elevation but not the 
crustal thickening could then be a consequence of thermal 
expansion of the underlying mantle. Although these various 
opinions diverge considerably, and the specific physical 
processes are vaguely described at best, the ultimate cause 
of Tibet's elevation, thick crust and volcanism is virtually 
unanimously attributed to the India-Eurasia collision. 

The existing data for Tibet are too scarce to discriminate 
among the various proposed evolutionary mechanisms, but 
tectonic models of the evolution of the India-Eurasia 
continental collision must. provide satisfactory explanations 
for the structure of the plateau. Furthermore, studying the 
geophys'^cal characteristics of the Tibetan plateau and the 
active Himalaya probably provides basic information that 
constrain the processes of large scale continental tectonics 
in general. 

The purpose of this study is to present a comprehensive 
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••Ismological data aat that placas In^rovad eonatralnta on 
tha aaiainic valocity atructura of tha crust and uppar nantla 
banaath tha Tlbatan plataau. Wa discuss saparataly tha 
four ralatlvaly Indapandant saismic nathods amployad: 
surfaca wava disparsion for paths crossing tha plataau# 
synthasis of a phasa crossing tha plataau# and 
rafraction profilas using aarthquakas in Tibat and tha 
Hiinalaya# and talasaismic S wava t^aval tima rasiduals from 
aarthquakas in and naar Tibat. In aach saction dascribing 
tha individual saismic phasas# a briaf introduction is givan# 
followad by tha rasults of obsarvations . Constraints on 
possibla models and discussions of tham ara than prasantad. 

At the end, we present an interpretation of all the 
observations and discuss their tectonic implications. 

None of tha studies of individual saismic phases requires 
a unique velocity structure, but a consistent irodol of the 
crustal thickness and uppermost mantle valocitias can be 
obtained from the combined observations. In particular, tha 
non-uniquenass of the velocity structure has bean reduced 
considerably bv the determination of tha uppermost mantle 
velocities from refraction profiles. 
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SURFACE WAVE DISPERSION 
a) Diaparalon Curvaa 

Wa analyzed long-period eeismograms for nine earthquakes 
in and near Tibet recorded at one or more of the four stations 
of the World Wide Stan^irdized Seismograph Network (WWSSN) 
just south of the Himalaya (Table I). More than 80% of each 
of the great circle paths is withiri the plateau. We chose 
not to use longer paths that cross much of Asia because of 
the difficulties in deciding how to regionalize them and 
because of the likelihood of marked lateral refraction at the 
boundaries of Tibet and elsewhere. We trust more the results 
obtained from paths confined much as possible to Tibet. 

Thus, the range of epicentral distances is limited to about 
20* (Pig. 1). The multi-filtering technique (Dziewonski et al., 
1969] vv'?.s used to obtain the group velocity dispersion curves. 

We analyzed only seismograms that have high signal- to-noise 
ratios and for which a minimum amount of multi-pathing was 
apparent (Fig. 2). For only two paths were clear long-period 
(^60 sec) Rayleigh waves recorded with large eunplitudes 
(Fig. 2a). Unfortunately Rayle.gh waves with periods near 
30 sec were so large on these two records that it was 
impossible to study the dispersion curves for periods less 
than 60 sec. 

We analyzed only fundamental modes of both Love and 
Rayleigh waves with periods from 10 to 90 sec. Group 


9 . 


velocity dispersion curves for Rayleigh waves differ from 
one another by about 3% (Fig. 3a) , which is comparable with 
that in some previous studies [e.g., Chun and Yoshii, 1977]. 
The dispersion curves between 30s and 60s are poorly 
defined because of a low signal-to-noise ratio, which is 
apparent from visual inspection of the seismogreuns. The 
amplitudes on the seismograms are reduced, at least in part, 
because the group velocity dispersion curve is very steep 
in this particular frequency band, but Bird and Toksdz 
[1977] inferred an extremely high attenuation in this band 
as well. Since the uncertainty of this part of the dispe ^ion 
curve is l^u:ge, we consider that this part is essentially 
undetermined. Compared with the previous studies mentioned 
above [e.g., Chun and Yoshii, 1977; Bird and Toksdz, 1977] 
our long-period Rayleigh wave group velocities are about 2 
to 3% faster at 70 to 80 sec but slower at longer periods by 
about the same amount. This difference contributes partly 
to the differences in preferred mantle velocity models. 

Love wave group velocities observed here also differ from 
one another by 5 to 10% (Fig. 3b) . No clear Love wave with 
period longer than 70 sec was observed. 

We also used *'hase velocities to aid in the discussion 
of these observations. Phase velocity dispersion curves 
(Fig. 4) are calculated using the standard single station 
method [e.g., Weidner and Aki, 1973]. 
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A combination of a mislocation of 10 km and a 2 iec 
error in origin time can cause an error of 0.8% in 
deteminations of both group and phaae velocities at a 
distance of 20*. An uncertainty in the group arrival time 
or phase delay of one-eighth of the period at periods near 
80 sec leads to an uncertainty of 1.6% in group or phase 
veloc'Lty at this distance. We assume a total combined 
maximum error of up to 3% in both is possible for the long 
period end of the dispersion curves. The uncertainty in 
the group velocity at short-periods should be less. However, 
lateral heterogeneity will have a much larger effect on 
relatively shorter wavelength components than on the longer 
wavelength components. This might be part of the reason why 
there is a fairly large scatter in the data at short periods. 

For shallow events, the source depth is not important 
at long periods. Because both events I and II have nearly 
vertical strike-slip fault plane solutions, obtained primarily 
from P wave first-motion polarities, [Molnar et al.. 1973; 
Tapponnier and Molnar, 19771 the uncertainty in the calculated 
phase delay is small. At extreme conditions, assuming an 
error of 0.2 cycles in various phase delays, an error of 3.4% 
in phase velocity is possible if all the errors conspired 
together for periods near 30 sec and a distance of 20%. 

The relatively large error in the phase delay assumed here 
takes into account the effect of a possible non-step-like 
source time function. 
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Not* that th« two Raylttlgh wav* aalamograms with largo 
•ignal-to-noita ratios at long psriods (1-NDI and II-SHL) 
Involvs different paths (Fig. 1) . The I-NDI path is almost 
entirely within Tibet except at the southwestern end where it 
crosses the Himalaya at a high angle. The II-SHL path crosses 
part of the Pamir before it enters Tibet and then meets the 
Himalaya at a low angle. The difference between the two 
observations at periods longer than 60 sec could be simply due 
to observational errors, but it might actually reflect real 
differences in average velocity structures for different patha. 
b) Velocity Models 

Resolving power of data and assumptions in models : 

Before interpreting the data, we briefly review the 
resolution of surface wave observations on velocity structures. 
Assuming a laterally homogeneous earth structure, Der et al., 
[1970] investigated effects of observational errors on the 
resolving power of surface waves using synthetic data with 
periods and distances comparable to those used here. 

They found that for similar period ranges the fundeunental 
Rayleigh mode yields considerably better resolution of the 
upper mantle velocity than the fundamental Love mode. 

If observations are restricted to fundamental modes only, 
crustal models with no more than 3 layers are justified. 

With both fundamental mode Rayleigh and Love waves and their 
first higher modes, the shear-wave velocity (Vg) resolution 
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in th« upper cruet, lower cruet, and uppennoet mantle are 
about 0.05, 0.08, and 0.1 km/sec, reepectlvely . In particular, 
the ehcar wave velocity reeolutlon la only 0.1 km/aec for a 
layer 100 km thick centered at a depth of 125 km. Unleaa the 
period range of group velocity obeervatlona for Rayleigh and 
Love wave fundamental modee exceede 80 sec, the upper mantle 
velocity reeolutlon la poor. The above summary Is appropriate 
for group velocities. Phase velocity observations provide 
somewhat leas resolution than group velocity observations with 
the same observational errors. 

Der et al. [1970] assumed standard deviations of group 
velocities of approximately 0.3% near 15 sec and 1.5% at 
80 sec; these errors are probably smaller than ours. With 
reliable group velocity observations of fundamental Rayleigh 
and Love waves at limited frequency band and with fairly 
large scatter (Fig. 3), and phase velocity measurements of 
fundamental Rayleigh mode between 60 to 90 sec only (Fig. 4) , 
we did not perform a formal Inversion. Instead, we exaunlned 
a range of plausible models; and then from a comparison with 
the observation, we assigned bounds to the range of possible 
structures . 

In constructing these models, we further assumed that: 

1. The velocity structure bolow 90 kra for most models 
and density structure for all models are the same as Gutenberg' 
continental model (.Oorman et al. , 1960] . 
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2. Tha ratio of P wave velocity to S wave velocity is 
fixed at 1.73 in the crust and 1.80 in the uppermost mantle. 

3. The shear wave velocity of the top 3.75 km is fixed 
at 2.55 km/sec [Chen and Molnar, 1975]. We did not consider 
the period range 5 to 10 sec here. 

4. No more than 3 layers in the crust* excluding the 
top sedimentary layer* are included and the mantle between the 
Moho and a depth of 90 km is assumed to be homogeneous. 

Since the resolution of P wave velocities is very poor 
compared with the resolution of S wave velocities* we examined 
dependencies of the dispersion curves only on V, in the mantle* 
the crustal thickness, and the shear wave velocity structure in 
the crust. Theoretical dispersion curves were compared 
primarily with Rayleigh wave group velocity observations* but 
for selected structures comparison was also made with Rayleigh 
wave phase velocity observations. Love wave group velocity 
observations served only as a loose constraint. As long as 
theoretical values fall within the scatter of these data* 
no further comparison with them was made. Surface wave 
dispersion at the short period range (<30 sec) does not 
constrain the crustal thickness or the upper mantle shear wave 
velocity much. Little attention will be paid to these short 
period range observations in the following discussion of 
these two subjects. 
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Crustal thickn««« and uppt mantl* S wav vlociti«s . 

Raylaigh wava dlaparalon curvas for thraa dlffarant 
crustal thlcknassas (SS, 70, and 85 km) wars calculatad for 
a ranga of upparmost inantla shaar wava valocitlas (0nt) from 
4.4 to 4.9 km/s (Figs. 5-7). At long periods, higher values 
of 8m causa highar Raylaigh wava group velocitias, and graater 
crustal thlcknassas cause lower group velocities (Figs. 5-7) . 
For a crustal thickness of 55 km, only the model with the 
lowest value of 8m km/s, model F4) fits the upper curve 

of the Rayleigh wave group velocity observations (Fig. 5) . 

A lower value of 8^ would be required to fit the lower curve, 
and if the assumed thickness of the crust is less than 55 km, 
8m must be lower than 4.4 km/s. For a 70 km thick crust 
(Fig. 6), however, the 8m seems to be constrained to be near 
4.7 km/s (Model S7) or higher in order to be consistent with 
observed group velocity dispersion curves. Values of 8m 
small as 4.4 km/s seem to be incompatible with a 70 km crust 
(Fig. 6) . If we further increase the crustal thickness to 
85 km (Fig. 7), the observations constrain 6m to be 4.8 km/s 
or even 4.9 km/s, (Fig. 7), but the fit between the observed 
and theoretical curves appears to us to be worse than for 
the previous cases of 55 km and 70 km thick crust. These 
results confirm Der et al.'s [1970] numerical experiments 
that show that surface wave dispersion alone is inadequate to 
constrain both the crustal thickness and the upper mantle 
shear wave velocity tightly. 
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Phase velocity dispersion curves for four velocity 
models are plotted with the two observed dispersion curves 
in Fig. 4. The curve for model F4 (55 km thick crust, 

■ 4.4 km/s) lies within the scatter of th'j observations. 
Both model S7 (70 km thick crust, ■ 4.7) and model SE 
(70 km thick crust, Sm " ^<8, with the lowermost crustal 
velocity lower than that of model S7) fit the observed phase 
velocities for the I-NDI path well and were within the 
uncertainties of the observations for the II-SHL path. 

Model S4 (70 km thick crust, > 4.4 km/s) does not fit 
the II-SHL observed phase velocities and is systematically 
lower than the data for I-NDI. It turns out that if 
8m 1 ^’8 km/s, however, the phase velocity data cannot rule 
out a crustal thickness of 85 km. 

Velocity structure of the crust; The crustal structure is 
constrained primarily by the shorter period portion of the 
dispersion curves. This portion of the dispersion curve 
depends primarily on the mean shear wave velocity in the 
crust and on its gradient. Fig. 8 shows the theoretical 
Rayleigh wave group velocity dispersion curves for 7 
different models with a fixed crustal thickness of 55 km 
and 6m • 4.4 )cm/s . If the acceptable values of the 
calculated Rayleigh wave group velocities are those that 
lie within the scatter of our data, then for such a crustal 
thickness, an average crustal shear wave velocity (8^) of 
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3.5 km/s to b« too high (Fig. 8 a) . A singl* layer 

batvraan the top low velocity surface layer and the Moho 
clearly will not fit the data unlesi s is much less than 
3.5 km/s (Fig. 8b). By adding a higher velocity layer at 
the base of the crust so as to cause a velocity gradient 
in the crust and maintaining a constant Q, dispersion curves 
can be constructed that agree well with the data (Fig. 8b) . 

To obtain such a fitr however, requires a relatively low 
velocity upper crust bo that 8^ is still less than 3.5 km/s. 

Calculated Rayleigh wave group velocity dispersion 
curves for 5 different crustal structures with a 70 km thick 
crust and 8,^ ■ 4.7 km/s, also are consistent with J being less 
than 3.5 )cm/s (Fig. 9). The long-period end of the dispersion 
curve however, is sensitive to the shear wave velocity in 
the lower crust. For instance, the addition of a relatively 
high velocity layer (such as Vg > 3.9 km/s) between 60 and 70 
km depth has the same effect on the long period portion of 
the dispersion curve that an increase of 8^^ (such as from 
4.7 to 4.8 km/sec) below 70 km has (Fig. 9b, c). These 
variations are, of course, only examples of many other 
possibilities. Thus there are trade-offs between the lower 
crustal shear wave velocity and 8,„ which ib in turn coupled 
with the crustal thickness. 
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For a crustal thicknass of 85 km, tha constraints on 
tha crustal structura placad by surfaca wavas will ba avan 
loosar chan tha pravious discussions for 55 and 70 km thick 
crust casas. Wa did not invaatigata any datailad roodals 
for this crustal thicknass. 

Tha most important constraint that surfaca wava disparsion 
places on tha crustal structura is that tha avaraga crustal 
shaar wava valocity is low. Claarly for tha uppar 40-50 km, 

's <3.5 km/s. Since the data seem to ba battar fit by a 

series of layers of increasing velocity than by a single 

layer of uniform velocity in tha crust below sadiroants, tha lower 

crust could have a shear wave valocity greater than 3.5 km/s, 

but it appears that even with this increase, is lass than 

3.5 km/s for the entire crust. Our data do not require a 

low velocity zone in the crust (cf. Chun and Yoshii, 1977]. 

The crustal thickness and uppermost mantle shaar wave velocity 
combinations obtained earlier are not very sensitive to the 
possible variations in the crustal structure, 
c) Synthesis of Rayleigh wave signals 

The group and phase velocity dispersion curves contain 
very little information about the amplitude of the surface 
waves other than the fact that there must exist a substantial 
amount of Fourier amplitude (spectral density) at a particular 
frequency in order to obtain a reliable determination of the 
group arrival time. As a further check on the velocity 
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■truetura# w« have usad tha auiranatlon of normal modai to 
ganarata ayninatlc aalamograma of vartical componanta of 
Raylaigh wavaa to Invaatigata tha poaslbility of datactlng 
structural dlffarancaa from thasa synthatlc aalamograma. 

Tha ganaratlon of aynthatic aalamograma raqulraa a knowladga 
of tha aourca dapth, tha fault plana solution, tha sourca 
tlma function, tha salsmlc momant, and paramatars dascrlblng 
attanuatlon and a valoclty structura. Usually, ona assumas 
a knowladga of savaral of tha factors and trlas to daduca 
tha othars fr«n tha agraamant of the synthatlc salsmogram 
with tha obsarvad one. The method was originally designed 
by Kanamorl (1970] to Investigate the long-period sourca 
characteristics of earthquakes at large distances. 

Wa suninad fundaunantal Rayleigh modos with periods between 
about 20 to 120s and generated synthatlc salsmogra.ns for 
I-NDI and ZI-SHL observations (Fig. 1 and 2) . Both tha 
synthatlc salsmogr 2 uns and tha observed seismograms (paths 
I-NDl and II-SHL, Fig. 2) were low-pass filtered with a cut-off 
frequency at 0.04 Hz (25 sac) to eliminate excess high 
frequency oscillations f'hat are superimposed on the long 
period signals of Interest. The finite filter length tends 
to Introduce different effects on the ends of the synthetic 
seismogreuns from those of the observations and Is a shortcoming 
for relatively short distances as in our cases. The abrupt 
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beginnings of the synthetic seismograms results from 
truncating the earlier part of the time series which is not 
of interest here. The following discussion deals with simple 
experiments of the effects of the different parameters on the 
time domain observations. 

Fault plane solutions for both events are held fixed 
and the uncertainties associated with them have small effects 
on the calculated wave form. A point source time history is 
assumed to be a step for simplicity. The severe attenuation 
for this region, proposed by Bird and Toksdz [1977] was applied 
to generate the theoretical seismogram marked with 'low Q* 
in Fig. 10a. The one marked 'rornal Q* was generated with 
the observed attenuation for Rayleigh waves in an average 
earth structure [summarized in Kovach, 1978]. Although the 
amplitude for the 'normal Q* synthetic seismograun is 351 
larger than that of the 'low Q' synthetic seismogreun, the 
wave forms are quite similar in the bandwidth of interest 
her'^, 60 to 90 sec (Fig. 10a). For both synthetic seismograms 
the same velocity model S7 (70 km thick crust, Sj^ ■ 4.7 k>/s, 
Fig. 6) and a source depth of 5 km were used. 

Although the observed wave form of the II-SHL path is 
also shown in Fig. 10a, the purpose of this figure is to 
illustrate the effects of a very low Q. Part of the mismatch 
between the synthetic and observed signals arises from the end 
effects due to a finite filter length. The coda of the S wave 
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and th« 20 mc p«rlod Raylaigh wavs ar« tha most noticsabla 
sourcss of nolss (Fig. 2a) . Ths calculatad Raylsigh wavas 
ara dalayad in tima with raapaet to tha obaarvad onaa in 
Fig. 10a. Tha maximvun mismatch in group arrival tima 
corraoponds to 2.7t diffaranca in group valocity, which agraas 
with tha disparsion curvas in Fig. 6 and is within tha 
obsarvational arrors. 

Fig. 10b showa thraa synthatic aaismograms, with thair 
paak-to-paak maximum anqplituda normalizad to a common siza. 
asauming modal S7 with diffarant aourca dapths of 5, 10 
and 30 km (tha actual amplitude ratio is about 3s2tl raspactivaly) . 
Tha diffaranca in wava forms for sourcas at 5 km and 10 km is 
small. Thara is a slight diffaranca in wava form for a 
sourca at 30 km dapth. This can ba attributad to tha 
intarfaranca of tha short pariod (<40 s) Rayleigh wavas, which 
hava a phasa shift of :r with raspact to tha longar pariods and 
much waakar amplitudas than for a dapth of 10 km (Tsai and 
Aki, 19701. 

Synthatic saismograms wara ganuratad for thraa diffarant 
valocity modals: ■ 4.4 km/s and 4.7 km/s for a crustal 

thicknass of 70 km (S4 and S7) and ■ 4.7 km/s for a crustal 
thicknass of 55 km (F7 ) (Fig. 10c). A 5 kJn sourca dapth and 
'avaraga-aarth* Q structure were used. With a constant 
Bgj (4.7 km/s), models with 55 km and 70 km thick crusts 
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(modal F7 and S7) gava qulta dlffarant aynthatic aaianograata 

(Pig. 10c) . Wa conaidar Modal P7 to ba Inadaquata to axplain 

tha obaarvation but that Modal S7 ia aatiafactory in matching 

tha obaarvad wava form of tha I-NDX path (Fig. 10c) . 

d) Summary of aurfaca wava raaulta 

» ~ 

1. From group valocity diapareion curvaa in tha pariod 
ranga of 10 to 30 aac. tha avaraga ahaar wava valocity in 
tha uppar cruat ia lowar than 3.5 )cm/a. 

2. From obaarvad and calculatad Raylaigh wava group and 
phaaa valocity diaparalon curvaa and from aynthaaizad 
Raylaigh wavaa with parioda batwaan about 60 to 90 aac, tha 
upparmoat mantla ahaar valocity cannot ba conatrainad wall 
aven with an aaaumad cruatal thic)cnasa. If wa aaauma a 8^ 

of 4.4 )cm/a, tha cruatal thic)cnaaa ia probably no mora than 
55 km. If ia about 4.7 km/a, tha cruatal thicknaaa ia 
likely to be greater than 65 km. A cruatal thicknaaa of 85 km 
ia not impoaaible if 6^ « 4.8 km/a and would be required if 
were yet larger. 

PHASE SYNTHESIS 

We observed one very clear Pjj phaac that croaaai Tibet 
(i'ig. 2a) . P^ propagatea aa a leaky mode in tha cruat and 
travels with velocity higher than the shear wave velocity of 
the underlying pu.ntle (Gilbert, 1964; Oliver and Major, 1961] . 
The P^ wave train arrives betw'<.«n the P and S phases and is 
normally dispersed. P^ can be viewed as the result of 
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multlpla reflected and rafractad P and SV wavas, aach of 
which la, at least once, an SV wave along the ray path 
[Ha: Jsargar, 1972]. 

Daapita the energy leakage, the form of the wave 
train can be used to infer the crustal structure. Both the 
phase and group velocity dispersion curves are fairly steep 
and the Pj^ wave train is generally short. This makes 
accurate measurements of dispersion difficult [e.g., Poupinet, 
1972]. Instead, we have synthesized seismograms of Pjj using 
different velocity structures in order to constrain the 
structure. Helmberger [1972] generated synthetic seismograms 
for Pj^ at distances between 4** to IS** from explosions using 
generalized ray theory, which employs a Cagniard-OeHoop 
algorithm. He assumed a layer over half-space. We followed 
his approach and used a double-couple earthquake source 
[Langston and Helmberger, 1975]. 

Vfe calculated synthetic seismograms oy summing 
generalized rays of up to suite order number 7 (7 segments 
in a ray) with a double-couple source at 5 km depth. The 
fault plane solution is the same as that used in the surface 
wave calculations before. No attenuation operator is 
applied to the synthetic seismograms. When we applied 
Futterman's [1962] constant Q operator with t/Q 1 sec 
[where t is the travel time) tlie difference in the waveforms 
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was baraly dlacernibla. (This corresponds to an average 
Q 300 for our case) . The assumed far-field source time 
history is given by t exp[-t/10 (sec) 1 » which would correspond 
to a fault length of about 30 km in Brune's [1970] formulation. 

At first, synthetic seismograms were generated for models 
with a layer over half-space, each of which represented an 
average of the velocities in models obtained from surface 
wave dispersion. Figures 11a and 11b show the synthetic 
seismograms for two models with 70 )cm and 85 )cm thick crusts 
respectively, together with the observed signal. The first 
large peak in the synthetic seisnKjgrams results from S to P 
and P to S reflections off the Moho. The large amplitudes at 
the beginning of the observed seismogram are the direct, P, 
pP and sP, which penetrated into the upper mantle to 400-600 km 
and were not synthesized here. Numerical experiments showed 
that variations in the mantle velocities between 8.2 to 8.5 
km/s have a negligible effect on the large amplitude pulse 
in the synthetic seismograms. Note that two quite different 
models (with a difference of 15 km in crustal thickness) gave 
similar wave forms in the synthetic seismograms (Fig. lib) . More- 
over the maximum difference in group velocity between the two 
synthetic seismograms is only about 2%. This suggests poor 
resolution of velocity structures with P^ observations alone. 
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N«xt> w« examined a eeries of multi-layered modela similar 
to those based on surface wave dispersion in the last section. 
Helmberger (1972] concluded that for a layer over half-space, 
the source depth is not crucial, as long as the source is in 
that layer. Also an additional top low velocity layer has a 
small effect on the long period character of the synthetic 
seismograms. From a series of numerical experiments with 
model SE (Fig. 9b) we found that the effect of source depth 
and that of a low velocity layer are insignificant as long 
as the source is in the layer next to the low velocity 
surface layer. We tried several models derived from the 
surface wave dispersion curves that included crustal reflectors. 
Only rays reflected between the free surface and the interfaces 
including and above the Moho were used. Fig. 11c shows one 
example of such a synthesis. This particular model SEMI has 
the same shear wave velocity structure as model SE (Fig. 9b) , 
except the upper crustal layer has a P wave velocity 7% less 
than that of model SE. The large Moho reflections seen in 
Fig. lib are suppressed. In addition, the late large arrivals, 
which resulted from the mid-crustal reflector, nave a 
significant effect on the long-period wave form (compare 
Figs. 11c and 11b). The synthetic seismogram from this 
three-layered cristal model fits the long period 
characteristics of the observation much better than the 
single layered crust cases. Other models were tested and 
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both volocltiss and layer thicknesses were modified in a 
trial and error process. No significant improvement was 
found with models of different crustal thicknesses and 
structures from that of model SE, and with a modified model 
S7 with a linear crustal-mantle transition. The model 
'TP-2', with a low velocity zone in the crust preferred by 
Chun and Yoshii [1977], does not provide a better fit either, 
even after some variations in its crustal and mantle structures. 
Presumably, with more layers in the velocity models we might 
achieve a better match to the observed P^ signal. Yet 
considering the poor resolving power on the structure, further 
tiral seems fruitless. It is likely that the gross form of 
P^ can be matched with crustal thickness of 55 to 85 km if 
there is a velocity gradient in the crust, but that the 
observations do not constrain the structure very much. 

REFRACTION PROFILES 

It is well-known in classical seismology that the P^^ 
and Sj^ velocities can be estimated from the travel time vs. 
distance plot using close earthquakes [e.g. Bullen, 1963, 
p. 194], This involves using many stations for each event. 
Although, to our knowledge, Lhasa (LHA) , has been the only 
permanent seismograph station in Tibet, the usual method 
can be modified using travel times from many earthquakes 
at a single station. 
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In practice, different earthquakes are reported with 
different focal depths. Moreover, the focal depths can 
easily be in error causing a corresponding error in the 
origin time. In either case, it is necessary to normalize 
the focal depth to a common depth in order to avoid introducing 
artificially large residuals due to different focal depths 
and therefore different path lengths. This can be dona by 
adjusting the origin time with a depth correction calculated 
from equation (3) in Bullen [1963, p. 194], assuming an 
approximate velocity structure. The difference between the 
assumed and actual origin time is assumed to be absorbed by 
this depth correction completely. Any incomplete trade-off 
between the focal depth and the origin time, together with 
a small error induced by an error in the assumed velocity 
structure will cause a possible systematic error in the 
measured travel time. Because this error is presumably 
Independent of distance in the normalized travel times, it 
should introduce only random errors in the evaluation of the 
velocity and therefore only affect the variance associated 
with the velocity determination, not the estimate of the 
velocity. 

and arrival times were taken from the Lhasa 
station reports available between 1960 and 1965 and again 
between 1971 and 1973 [Institute of Geophysics, 1966, 1974]. 
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We eelected well- located events, for which nK>re than 35 
readings from teleseismic stations were used in the location. 

We used the epicenters, origin times, and depths given in 
the Bulletins of the International Seismological Summary (ISS) 
and the International Seismological Center (ISC) . 

across Tibet; All 61 reported arrival times between 3* 
to 17* from Lhasa are plotted in Fig. 12. A model with a 
crustal layer with P wave velocity (Vp) of 6.1 km/s over a 
half-space with Vp ■ 8.0 km/s was assumed to normalize the 
focal depth to 33 km, which is near the mean reported depth. 

For all cases, we took the radius of the interface to be 
6306.2 1cm. For a Moho at 35 km depth, instead of 70 km, 
the estimated velocity should be increased by 0.5%. The P^ 
velocity is determined to be 8.15 + 0.04 km/s, compared with 
8.27 + 0.07 km/s if no focal depth correction is applied 
(Table II). If the apparent velocity were in fact 7.9 )cm/sec, 
then the data should be more nearly parallel to the dashed 
line in Fig. 12. Such a low Pj^ velocity will not fit the data 
well unless there is some special distance dependent systematic 
error in locations, depths c*“ origin times. If we limit 
ourselves to sharp Pj^ arrivals (excluding those reported as 
eP or ePn in the station reports) , we obtain a P^ velocity 
of S.20 + 0.04 km/s from 42 arrivals. Without the focal 
depth correction, we obtain 8.24 + 0.09 km/s (Table II). 
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If w« consider only the 39 sharp arrivals between 3* 
to 10* (Fig. 13), the P„ velocity Is estimated to be 8.12 * 

0.06 km/s for the depth corrected case (Fig. 13b). This Is 
the only case In all the profiles discussed In this section 
for which the velocity determined without a depth correction, 
8.0 * 0.1 km/s (Fig. 13a) Is thf; smaller. A total of 19 
emergent arrivals between 3* to 17* gave an estimate of the 
velocity of 8.18 * 0.13 km/s. Thus the entire data set is 
very consistent (Table II) . 

Pn across the Himalaya: The data discussed above Includes 

only events north of the Indus Tsang-Po suture because of the 
possibility of a different structure beneath the Himalaya. 
However, recordings at Lhasa for 24 events south of the suture, 
for nearly all of which sharp readings were reported, yields 
a Pn velocity of 8.11 > 0.07 km/s. The epicentral distances 
range from 3* to 17*. Among these, 21 events occurred within 
10* of Lhasa (Fig. 14). Despite the scatter, the same estimate 
for Pn velocity is obtained. Thus no significant difference in 
the Pn velocity is observed for regions north and south of the 
suture zone. 

Sn across the Tibet and the Himalaya: There are also 23 sharp 

Sn arrivals reported at Lhasa from events north of the suture 
in the distance range 3" to 10®. (Only three such arrivals 
were reported tor distances greater than 10®). After the depth 
normalization, the Sn velocity is estimated to be 4.8 + 0.1 km/s 
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(Fig. 15). There is a much larger scatter in the S„ residuals 
than the residuals probably because of the slower 
velocity and because of difficulties in identifying the 
arrival on the seismograms. A low velocity of 4.5 km/s 
would not fit the data well (Pig. 15). The 11 sharp arrivals 
for events south of the suture, reported at Lhasa within 3* 
to 10* gave the same velocity (Table II) . 

Intercept times and crustal thicknesses; Despite the large 
uncertainties, the intercept times for all suggest the 
depth of the interface (the Moho) to be 65 + 30 km. 
intercept times have even larger uncertainties, but the depth 
of the interface calculated from them fall within the range 
(Table II) of 80 ^ 50 )cm. The epicenters of all the events 
used in this refraction study show a broad azimuthal coverage 
and distribution of epicentral distances with respect to Lhasa 
(Fig. 16) . Nevertheless, because purely reversed profile 
cannot be made, it is remotely conceivable that a dipping Moho 
or a systematic horizontal variation in the velocity structure 
centered around Lhasa might render the real mantle velocities 
different from the measured apparent velocity. For models 
with a crustal layer with Vp ranging from 5.9 to 6.5 I«m/s 
overlying a mantle half-space with Vp ranging from 7.9 to 8.4 
Icm/s, a uniform dip of about 0.8® of the Moho could add an 
uncertainty of +0.1 km/s to the estimated velocity. 

Such a dip would correspond to a difference in crustal 
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thlckn«st«a of 15 km ov«r 10* dlatanc*. Tharafora, if Lhasa 
wara to ovarlla a maximum dapth to tha Moho baneath Tibat, 
such that tha crustal thicknass dacraasad by 15 km, 10* to 
both tha north and wast, than tha valocity could ba 
ovarastimatad by 0.1 km/s. Wa think that bacausa tha avaraga 
alevation of Tibat is so uniform, a drastic variation in 
crustal thicknass is unlikaly, and it would ba fortuitous 
if Lhasa wars to overlia tha daapast crust or tha cantar of 
a systematic horizontal variation of tha velocity structure. 
Summary of refraction results; High values of tha P^j and Sj^ 
velocity beneath Tibet are obtained from refraction profiles 
using earthquakes in Tibet and arrival times at Lhasa. 

Assuming wi Moho at 70 km depth, the P^^ velocity is estimated 
to be 8.12 km/s with a minimvun formal standard deviation of 
0.06 km/s covering the region north of the suture from 3* to 
10* from Lhasa. No significant difference in Pn velocity was 
found for the regions north and south of the suture. With a 
somewhat larger uncertainty, the velocity of is similarly 
determined to be 4.8 + 0.1 km/s. Both P„ and velocities 
are about the rame as those of the shields and stable platforms. 
All the intercept times have large uncertainties and they 
indicate a depth of the Moho of 65 -f 30 km. 
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TELESEISMIC S>P TRAVEL TIME RESIDUALS 

To our knowl«dg«, the high velocities for the uppermost 
mantle under Tibet discussed in the last section have not 
been proposed before. The vertical extent of the high 
velocity zone is of great importance in the interpretation 
of the tectonics. The refraction profiles require a layer 
beneath the Moho with a thickness of only about 10 km. Pjj 
is relatively insensitive to mantle properties [Poupinet, 1972] 
and the surface wave dispersion curves are not very sensitive 
to the velocity structure below about 100 km. The teleseismic 
travel time measurements discussed in this section are intended 
to serve as a crude estimate of the vertically averaged velocity 
structure to a greater depth ('^300 km) . 

Teleseismic P and S wave arrival times were measured from 
the records of the WWSSN stations for 9 earthquakes in Tibet 
and 8 earthquakes in the Himalaya (Fig. 17, Table III). Arrival 
times were measured only for S wave signals that are clear on 
the long-period records (Fig. 18). Uncertainties are about 
one second. To avoid inclusion of converted phases (e.g., Sp) , 
we tried to measure only the transverse components of the S 
waves. S-wave polarization angles were checked with known P 
wave fault plane solutions (Molnar and Tapponnier, 1978; 

Molnar et al., 1973, 1977; Ni and York, 1978]. We considered 
only the distance range between 30® and 80® to avoid 
triplication in the S wave travel time curve, and possible 
misidentif ication of SKS for S. 
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Both P and S wava traval tima realduala, with raapact 
to tha Jaffraya-Bullan [1940] (J-B) traval tlma tablaa* 
wara plottad aa a function of apicantral diatanca uaing 
locationa and origin timaa givan by tha XSC. Sinea tha 
hypoeantara and origin tiitiaa raported by ISC wara calculatad 
primarily uaing tha J-B P wava traval timaa » tha P wava 
raaiduala wara automatically minimizad. Bacauaa lata (or 
aarly) S wavaa ara likaly to ba aaaociatad with lata (or 
aarly) P wavaa, tha S-P raaiduala maaaurad hara probably 
giva lower bounda to tha S wava raaiduala. 

For 9 earthquakaa in Tibat, P wava raaiduala acattar 
around tha maan aa expected from tha hypocanter location 
procedure, whereas S wava travel times show large delays 
(Fig. 19). The S-P travel time residuals have correspondingly 
large delays (Table III, Figs. 19 and 21). Applying station 
corrections (Sengupta, 1975; Sengupta and Julian, 1976] to 
these data does not significantly reduce the scatter, but 
increases the S-P travel time residuals by about 1 sec 
(Table III, Figs. 19 and 21). 

With the seune procedure, we studied 8 earthquakes in 
the Himalaya (Fig. 20) and observed smaller S-P travel time 
residuals than for earthquakes in Tibet (Table III, Figs. 

20 and 21) . The rays from earthquakes in Tibet and the 
Himalaya to a given station more than 30** from the epicenters 
will have nearly identical paths except near the source 
region. Thus we infer that the large S-P travel time 
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d«l«ys for oarthquakos In Tibat ara dua to tha madlum 
propartias of tha ragion banaath Tibat and not to ayatamatlc 
arrors in tha traval tima tablao or to tha madium along tha 
raat of tha ray paths. 

In ganaral. thara is a strong trada-off batwaan tha 
focal dapth and tha origin tima for avants locatad 
talasaismically . Erronaous dapths could also, at laast 
partially, ba raaponsibla for tha larga S-P traval tinM 
dalsys for aarthquakes in Tibet, if thasa aarthquakas actually 
occurrad at a shallowar dapth than was obtainad in tha location 
procadura. Tha S-P traval tima rasiduals wara calculatad for 
each event assuming three possible dapths: their dapths given 
by the ISC, 5 km, and 33 km (Fig. 21) . To make an axtrama 
comparison, assume all earthquakes in Tibat occurrad at 5 km 
dapth, but all earthquakes in the Himalaya occurrad at 33 km 
depth. On the average, there is still a difference of 
approximately 2 sec between the S-P residuals from earthquakes 
in these two regions (Fig. 21, Tedsle III). 

Thus the average S-P travel times are late for avants 
in Tibet not only with respect to the J-B tables which 
represent an averaged earth model but also with respect to 
those for events in the Himalaya. Because low P wave 
velocities probably cannot be detected through routine 
location procedures using the P wave arrival times alone, 
a delay of Is in P wave travel times from the average earth 
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could oaBily axist. This would laad to anothar sacond in 
daisy in tha 8 wsva traval tima. 

Tha maaaurad S-P traval tima rasiduals. roundad to 
tha smallar naarast intagar. ara givan in 1 sac intervals 
in Figure 17 together with tha locations of tha 17 earthquakes. 
For this figure, tha focal depths for Tibetan events ara 
normalized to 5 km and for Himalayan events, at 33 km. 

Although our data points certainly do not allow a meaningful 
contouring of these residuals in Tibet, tha two earthquakes 
in tha canter of Tibet have especially large S-P traval time 
residuals (>5 sec), and the distribution of S-P travel time 
residuals in the Himalaya is relatively uniform. The 
immediate inference from the late S-P travel times in Tibet 
is that the relatively high P^ and Sj^ velocities obtained from 
tha refraction profiles discussed in the last section are 
not likely to extend to a great depth beneath Tibet. 

SEISMIC WAVE VELOCITY STRUCTURE 

Observed group and phase velocity dispersion curves 
of Rayleigh waves provide constraints on crustal thickness 
and the uppermost mantle shear wave velocity. For possible 
thicknesses of 55 km to 70 km, or even 85 km, the corresponding 
uppermost mantle shear wave velocities of 4.4 km/s nt less 
to 4.7 'V 4.8 km/s, or even as large as 4.9 km/s, are required 
respectively. We consider an velocity of 4.9 km/sec to 
be unlikely and therefore conclude that the crust probably 
is not an thick as 35 km. 
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Th« paths for tha obsarvatlons of surfaca wavas and tha 
phaca al.\ cross tha Himalayas. This apparant lataral 
inhomoganaity introducas an additional uncartainty that is 
difficult to quantify. Tha two claar long-pariod Raylaigh 
wava signals studiad hara cross tha Himalaya at vary diffarant 
anglas (Fig. 1) , and tha 3% diffaranca in tha group valocitias 
could ba a consaquanca of lataral hataroganaity . Nonathalass# 
wa cannot confidantly attributa this diffaranca solaly to 
raal structural diffarancas becausa of tha uncartaintias in 
the group valocity detarminations. 

Sinca both tha and valocitias wars astimatad to ba 
ralativaly high from tha rafraction profilasi Pn valocity 
8.12 + 0.06 km/s, velocity 4.8 + 0.1 km/s, tha crustal 
thickness is apparently about 70 km. It is difficult to 
attach a formal uncertainty to this thickness i but from the 
comparison between the observed and theoretical dispersion 
curves, we estimate that the crustal thickness is between 
65 km and 80 km (Figs. 5, 6, and 7). These models seem to 
be compatible with the Pl wave form, but it does not provide 
a strong constraint on the velocity structures by itself. 

Travel time measurements for earthquakes in Tibet to 
utations at distances of 30** to 80° indicate late S waves 
and long delays between P and S. Since the Tibetan earthquakes 
have late S-P travel times compared with earthquakes in the 
Himalaya, the average velocity beneath Tibet must be lower 
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than that of tha Himalaya. The relatively high velocity zone 
beneath the Moho of Tibet, estimated from the refraction 
profiles therefore probably is not thick. 

The S-P travel time residuals can be used to estimate 
differences in average upper mantle velocities beneath Tibet 
and stedale platforms or shields with some additional assumptions. 
Suppose all such differences were confined to the upper 240 km 
as inferred for shields and ocean basins [Dziewonski, 1971; 

Okal, 1977]. Also assume that the measured S-P travel time 
residuals are estimates of S wave travel time residuals. Let 
us compare the Tibetan models with model S2 of Dziewonski 
[1971] for the shield regions, but modified to include 70 km 
thick crust. As an extreme let us assume that the earthquakes 
in Tibet occurred at a depth of only 5 km. Then if one 
increases the velocities of model S2 (with a 70 km thick crust) 
by about 1% to correct for the effect of velocity dispersion 
due to attenuation [Hart et al., 1977], the vertical one way 
S wave travel time between the depths of 0 and 240 km for 
modified model S2 is 53.8 sec. This is about 1.1 sec faster 
than that for the Jef freys-Bullen model. Therefore the 
observed average S-P travel time delay of 3.6 sec with respect 
to the J-B tables for earthquakes in Tibet corresponds to an 
average difference of 4.7 sec in S wave travel time from the 
modified model S2. Of this 4.7 sec, 2 sec can be attributed 
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to tile offoct of a thick crust. The resulting S wav* travel 
time difference of 2.7 sec corresponds to a reduction of the 
shear wave velocity of 6% between 70 to 240 kro from that of 
model S2. If we assume that tnodels S7, SE (70 km thick crust. 
^ ■ 4.7 -4.8 km/s, Figs. 6 and 9) or F4 (55 km thick crust, 
fijn • 4.4 km/s. Fig. 5) were appropriate for Tibet in the 
depth range from 0 to 100 km instead of S2. then the average 
shear wave velocity beneath Tibet must still be about 4 to 5% 
lower than that of model S2 tor shields between the depths 
of 100 km and 240 km. Insofar as Dziewonski's [1971] model 
52 is appropriate for the shields, this would correspond to 
the average shear wave velocity of about 4.3 km/s between 100 
and 240 km depth beneath Tibet. Because of S-P residuals 
are lower bounds to the S wave delays and because we assumed 
extreme for the depths of the events in Tibet, the average 
shear wave velocity beneath Tibet could be still lower. 

The effect of velocity dispersion due to attenuation 
has been neglected in the above discussions of crustal 
thickness. It is difficult to estimate the contribution 
from such velocity dispersion without sufficient data on 
both attenuation and phase velocity [e.g.. Lee and Solomon. 
1978]. Qualitatively, if strong attenuation exists, the 
estimate of shear wave velocities deduced from surface 
waves should be increased by a small amount ("^1%) above 
those considered here before comparing them with body wave 
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observations [Hart et al., 1977]. Thus thera might ba a 
small systamatic ovarastimate of a faw km of tha crustal 
thicknass, but this ovarastimate is probably much lass than 
tha uncertainty quoted above. 

TECTONIC IMPLICATIONS 
A. The Thick Crust 

The apparently very thick crust ('v70 km) beneath Tibet 
implies that the crust has been thickened approximately by 
a factor of 2 since the collision between India and Eurasia. 
This is consistent with the idea of horizontal compression 
and shortening (Dewey and Burke, 1973]. The hypothesis of 
underthrusting the Indian plate beneath the Eurasian plate 
can be argued as still possible if somehow either the crust 
of Tibet were detached from its underlying mantle while the 
Indian plate was subducted along the base of the crust or 
if the Indian crust plunged into the mantle and then migrated, 
probably as a melt, up to the base of Tibet's crust. In 
other words, continental crust of the order of 100 to 1000 km 
length would have been consumed and recycled back near the 
surface in a time span of about 40 m.y. The current seismic 
data cannot discriminate between these two extremes but the 
relatively low avera'7e velocity of the crust implies that the 
crust is relatively hot. If cold Indian lithosphere were 
thrust under Tibet, it would have to have been heated since 


then. 
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B. Uppermost Mantle Seismic Velocities 

Both and velocities indicate a high velocity 
uppermost mantle under Tibet. The late S-P travel time 
residuals for earthquakes in Tibet indicate that these high 
velocities do not reflect the properties of the deeper part 
of the upper mantle. Although the interpretation of the 
travel time residuals are higlily non-unique, this high 
velocity layer extending no deeper than about 100 km under 
Tibet is consistent with the observations. 

Implications for the temperature in the mantle; The 
apparently high uppermost mantle velocities under Tibet are 
not directly comparable to those under the stable platforms 
and shields because the extremely thick crust in Tibet will 
create about 10 to 12 kb of additional lithostatic pressure 
at the base of the thickened Tibetan crust. Assuming that 
the pressure and temperature dependence of seismic wave 
velocities for olivine are representative of the mantle 
(e.g. Anderson et al . , 1968; Birch, 1969], 10 kb of pressure 
causes an increase in Vp of about 0.1 km/s and a temperature 
increase of about 250‘*C causes a similar decrease. has 

very similar pressure-temperature dependence to that of the 
Vp for most olivine group minerals. Exceptions are the 
forsterite samples measured by Anderson et al. [1968] , which 
have been questioned by Birch [1969]. Thus, if the and 
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velocities were tne same beneath stable platforms and Tibet 
and if the temperature and pressure dependencies of the 
seismic velocities were those appropriate for olivine, 
then the temperature of the Moho beneath Tibet would be 
about 250 *C hotter than that of the platform. This estimate 
of relative temperature would have an uncertainty as large 
as 120'*C, even if the resolution of the difference in 
velocity for Tibet and shields were only +0.05 km/sec. 

Thus this suggested difference in temperature is only a 
crude estimate, not a precise measurement. 

The temperature under platforms and shields are difficult 
to estimate. From heat flow measurements, Sclater et al. 

[1980] estimated a lower bound on the temperature for the 
Canadian shield to be less than 400 "C at the base of the 
crust. On the other hand, if we take the measured reduced 
heat flow as an upper bound of the mantle heat flow for 
Eastern U.S. (0.80+0.05 HFU, (Lachenbruch and Sass, 1976]) 
and the Indian shield (0.93 HFU, [Rao et al., 1976)), the 
mantle heat flow alone gives 470'’C and 540°C at the base of 
the crust respectively, without considering the effect of 
radioactivity and assuming a crustal thickness of 35 )cm and 
thermal conductivity of 6 x 10 ^ cal/cm®C-sec . Uncertainties 
in conductivity and heat flow measurements alone will probably 
add 20% uncertainties to those estimates. In any case, if 
the temperature at the base of the crust beneath platforms 
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and shialds were about 500 *C, then at the Moho beneath Tibet 
it could be about 750”C. Although the uncertainty in this 
number is very large, our point is that the high and 
velocities do not require an unusually cold upper mantle. 
Implications for the strength of the mantle ; At first glance, 
the relatively high and velocities are likely to be 
taken as evidence for a cold and therefore strong upper 
mantle. Because crustal shortening requires strong deformation 
of the underlying mantle, the inference of high strength 
might be taken as evidence contradictory to crustal shortening. 
In this section we contend that the current uppermost mantle 
is not very strong. 

The temperature at the Moho strongly affects on the 
mechanical behavior of the lower crust and upper mantle. 

One can estimate the mechanical strength of the uppermost 
mantle beneath Tibet assuming that the flow laws for steady 
state creep of olivine are representative of the uppermost 
mantle rheology [Goetze and Evans, 1979; Kirby, 1977; 

Tapponnier and Francheteau, 1978]. For differential stress 
(01-03) less than 2 kb, the flow law can be written in the 
form: 


,01-03 


''e.xp(-Q/RT) 


€ » 


( 1 ) 
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wh«r« *e is the strain rate, A and n are constants, Q it the 
activation energy, R is the gas constant and T is the absolute 
temperature. For dislocation creep of olivine, Goetze [1973] 
reported n ■ 3, A ■ 70 (bar“^ sec”^) , and Q ■ 122 (kcal/mole) . 
Note that the uncertainties in the estimate of activation 
energy and temperature will dominate the uncertainty of the 
estimate of the differential stress or strain rate. Assuming 
a strain rate corresponding to 100% strain for the past 40 m.y., 
since the collision, the calculated differential stress is 
1.1 kb at 750“C or 430 bars at 800“C. If the temperature were 
as low as 700**C, the calculated differential stress would be 
larger than 2 kb. Post [1977] obtained a quite different 
olivine flow laws for Mt. Burnet dunite. Using his 'dry' 
flow low under the same assumptions, differential stress is 
estimated to be 3.7 kbars at 800“C while the 'wet' flew law 
requires only 30 bars of differential stress under these 
conditions. If Coble creep (syntectonic recrystallization), 
instead of dislocation creep, however, is the predominant 
deformation mechanism, the calculated differential stress is 
only 210 bars at 800®C and 430 bars at 700®C at the given 
strain rate (Goetze, 1978]. Thus the lower bound of the 
differential stress required to deform the uppermost mantle 
under Tibet at the given strain rate through steady state 
ductile flow of olivine is only about 0.4 kb. Although the 
estimate of temperature beneath Tibet is too uncertain to 
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prove that the uppermost mantle is not strong, the apparently 
high and velocities beneath Tibet do not require high 
mechanical strength for the uppermost mantle. 

Implications on volcanism in Tibet; It is interesting to 
note that at temperatures of about 800^C and 20 kb pressure, 
incipient melting of crustal material (but not the uppermost 
mantle) is possible in the presence of a small amount of water. 
Fig. 22 shows solidi for some likely crustal and upper mantle 
materials based on data from experimental petrology [Lambert 
and Wyllie, 19^'0; Merrill et al., 1970; Merrill and Wyllie, 
1975]. None of the data have error bars and therefore they 
can only be used qualitatively. Nevertheless, the solidus for 
water saturated peridotite is over 1000**C despite the pressure 
(Merrill et al., 1970]. With only 0.1 weight% of water, 
partial melting of gabbroic material is possible at 20-25 kb 
pressure and between 700 to 800®C [Lambert and Wyllie, 1970]. 
The solidus for peridotite with 0.1 weight! water at this 
pressure range is over 1250“C. Thus for Tibet, there could 
be a geotherm such that a partially molten crust overlies a 
subsolidus uppermost mantle. The inference of a high velocity 
mantle lid on the low velocity zone, then, is not incompatible 
with the widespread volcanism over Tibet unless the volcanic 
activity is predominantly of uppermost mantle origin or unless 
the lower crust of Tibet consists of only anhydrous minerals. 
The geochemical evidence gathered so far is not sufficient 
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to roach any conclusion on whether the lower crust of Tibet 
is hydrous or not. Simply based on the fact that no 
tholeiitic basalt has been reported [Hennig, 1915, Deng, 

1978], it is believed that the volcanic activity involved 
mostly melting of crustal materials (e.g. Dewey and Burke, 

1973]. However, the unusually low average crustal velocity 
for Tibet is also consistent with high temperatures and melting. 
C. A Discussion of the Thermal Evolution of the Plateau 

The presence of widespread volcanic activity and diffuse 
seismicity with prominent normal faulting in Tibet (Molnar 
and Tapponnler, 1978) all indicate that Tibet is tectonically 
active at present. Yet the geothermal gradient in the crust 
shtuld have been decreased by crustal thickening. Thus the 
source of heat deserves some consideration. There are three 
apparent heat sources that contribute to Tibet's crustal 
geotherm: strain heating from the thickening of the crust, 
radioactive heat production in the crust and the mantle heat 
flux from below the crust. 

Since the crust of Tibet seems to be twice as thick as 
normal continental crust, the average strain rate in the 
crust is likely to be of the order of 1/40 m.y. Assuming 
an averaged differential stress of 1 kb throughout Tibet, 
the stress heating rate per unit volume (stress x strain rate) 
is only about 8 x lO"^ erg/cm^-sec or 0.2 x lO”^^ cal/cm^-sec 
(0.2 HGU) . This is considerably smaller than the presumed 
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radioactive heat production rate of about 1 HGU for lower 
crustal materials [e.g. Tilling et al. , 1970] and can be 
neglected. 

We investigate the effects of the mantle heat flux and 
radioactive heating in the crust by considering the simple 
one-dimensional heat conduction problem with the following 
assumptions : 

1. A constant mantle heat flux (Q^) is held at the base 
of the crust. 

2. The thickening of the crust is instantaneous emd 
the original geotherm is stretched accordingly. This serves 
as the initial condition. 

3. Tibet defines <i heat flow province such that the 
empirical relationship 


Qs - Qh * D Aq, (2) 

between surface heat flow (Qg) and surface heat production 
rate (A^) is valid. 

4. There has been no erosion. This assumption is 
justified by the facts that the Cretaceous limestone is still 
widesf''read in southern Tibet (Hennig, 1915, Morin, 1946], 
and that the present-day drainage system in Tibet is internal. 
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W« teparate tha tempera tur a T in the crust into the 
contributions from tha mantle heat flux and from 
radioactive heat production Tj^. Three possible distributions 
of radioactive heat producing elements in the crust are 
considered separately: (a) a uniform distribution with depth, 

(b) an exponentially decaying distribution as a function of 
depth with the surface heat production rate unchanged by 
crustal thickening, i.e., with no migration of heat producing 
elements after crustal thickening, and (c) an exponentially 
decaying distribution with depth but with the thickness of 
the layer of the variable heat producing elements unchanged 
by crustal thickening (i.e., the migration of radioactive 
heat producing elements is much faster than heat conduction) . 

The detailed derivations are given in Chen [1979]. We will 
only summarize the most important results here. Throughout 
this discussion, the crustal thickness is taken to be 
thickened by a factor of two to a value of 70 km. We also 
assumed that the thermal conductivity (k) « 6 x 10"^ cal/cm-®C-sec 
and the thermal diffusivity (<) « 0.01 cm^/sec. , and we shall 
concentrate our discussions only on the temperature at the 
base of the crust where it will be highest in the crust. 

First, consider only Tj^ with two starting temperatures 
(Tjn®) 500“C and 600*C, corresponding to mantle heat flow {Q^) 
of 0.86 cal/cm^-sec (HFU) and 1.03 HFU, respectively, at 
steady state. will reach 800*C in about 45 and 14 m.y., 
respectively. 
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Thtt contribution from radioactiv* boat production 
strongly dapanda on tha distribution of thosa alamants, 
paramatarizad by tha constant 0, in aquation (2). Tha 
staady stata tamparatura dua to radioactiva hast production 

is approximataly proportional to D^. For an axponantial 
distribution of radioactivity « D must ba much smallar than 
tha crustal thicknass to satisfy aquation (2) [Lachanbruch> 
1970]. If the differentiation of tha radioactive heat 
producing alamants is mucn faster than heat conduction so 
that D does not change during crustal thickening, than Tf^ 
contributes little to tha crustal gaotharm. This occurs 
because in most models, and probably in the earth also, 
the source of heat far from tha base of tha crust. 

For cases (a) md (b) , Tj^* k, 4T^* and after 40 m.y., 

^ contribution to tha initial 

geotherm due to radioactivity. For instance, at 70 km depth, 
for case a, with ■ 0.73 HFU, D ■ 25 km and A(z) ■ Aq ■ 1 HGU 
(T® ■ 490®C) , in 40 m.y. T ■ 801®C. The steady state final 
heat flow at the surface would be only 1.25 HFU for this case. 
Since heat flow measurements are not available for Tibet, 
these hypothetical values are only meaningful as plausibility 
arguments, and they should not be taken as fact. Nonetheless, 
the above calculations suggest that it is possible that the 
lower crust of Tibet could reach 700-800®C in 40 m.y., 
simply from the recovery of the geotherm provided that the 
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h«at flux from th« mantl* la graatar than about 0.9 HFU 
ragardlaaa of tha haat production in tha crust. Tha starting 
tamparatura at tha basa of tha crust is crucial in datarmining 
whathar tha tamparatura in tha lowar crust raachas 700 to 
800*C in 40 m.y. With a raasonably warm starting tamparatura 
of 490*C at tha Moho. and if thara is sufficiant radioganic 
haat production in tha crust to produca about 0.5 HFU aftar 
tha crustal thickaning/ then, for a mantle heat flux of only 
0.75 HFU, tha calculated temperature at the Moho could reach 
801*C aftar 40 m.y. Thus radioactive haat production could 
be important for the recovery of crustal gaotharm in 40 m.y., 
especially if radioactive heating also makes a significant 
contribution to the original geotherm [Bird et al., 1975]. 
Similar conclusions on the importance of ratioactive heating 
are reached by England [1979] for the thickened Archaean crust. 

Crustal shortening beneath Tibet by horizontal shortening 
[Dewey and Burke, 1973] seems to be compatible with both 
the seismic observations and these thermal calculations. Yet 
the behavior of the lithosphere as a whole during the 
continental collision is not clear. In fact, the assumptions 
in these calculations are directly related to this question. 

If we assume that at 70 km depth the mantle heat flux is 
about 0.9 HFU and that the temperature at the Moho is about 
800®C, then for a coefficient of thermal conductivity (k) 
of 7.5 X 10“^ cal/ cm/sec*C in the mantle, the 1300**C isotherm 
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will b« rtachcd *t about 112 km. If indaad tha lithoaphara 
wara maintalnad at an approxlmataly conatant thlcknaaa by a 
haat flux of about 0.9 KPU at tha baaa of tha non- convac ting 
layar* than tha aatimatad tamparatura at tha baaa of tha 
Tibatan cruat could ba axplainad aimply in tarma of tha 
racovary of tha gaotharm with littla contribution by 
radioactivity. Tha volcanic activity, at laaat in part, 
could ba axplainad by tha dapraaaion of tha lowar cruat to 
depths whara tha tamparatura ia high anough for it to malt. 

No anomaloua haat flow from tha as thenosphara would ba raquirad. 

Prasumably small seals convaction maintains an assantially 
constant haat flux at tha base of tha plats. Parsons and 
McKenzie (1978] describe this in terms of tha stability of a 
thermal boundary layer above tha viscous asthenosphara and with 
a thickness 6 between two isotherms. If a critical Rayleigh 
number Ra^ “ is exceeded, tha boundary layer becomes 

unstable so that small scale convection occurs and provides 
convective heat transport. In Ra^, g is the gravitational 
acceleration, a the coefficient of volumetric thermal expansion 
and V is the kinematic viscosity. If the whole lithosphere 
with the thermal boundary layer were thickened by horizontal 
shortening, 6 would increase and Ra^ would increase as 5'. 
Therefore crustal shortening and thickening would be a very 
effective way of enhancing the instability of the boundary 
layer and causin small scale convection. Effectively there 
would be a detachment of the bottom of the thickening 
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lithoaph«re ^hat might tend to keep the lithosphere from 
continuing to thicken while the crust does. 

CONCLUSIONS 

1. The crust beneath Tibet is about 65 to 80 km thick 
with an average shear wave velocity of less than 3.3 km/s 
from surface wave dispersion/ refraction profiles and 
Pj^-phase synthesis. 

2. Pjx and S^ velocities were estimated to be 8.12jk0.06 
km/sec and 4. 8+0.1 km/sec respectively. 

3. The averaged vertical shear wave velocity between 
70 to 240 km depth is about 6% lower than that of the shield 
model S2 [Dziewonski, 1970] and this average velocity beneath 
Tibet is also lower than that below the Himalaya by about 
the same amount. 

4. From the and velocities and the pressure- 
temperature dependence of the seismic velocities for olivine, 
the temperatures at the base of the Tibetan crust could be 
about 250*’C hotter than that beneath shields or stable 
platforms. Since the temperature at the Moho beneath platforms 
is about 500 *C, the base of the crust beneath Tibet could he 
near 750‘*-800"C. 

5. At this temperature, and depending on the deformation 
mechanism, the uppermost mantle is likely to have a maximum 
mechanical strength of less than about 1 kb. This value 
decreases rapidly with depth in the uppermost mantle. 
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6. Partial melting of the lower crustal materials but 
not the uppermost mantle is likely with these estimated 
temperatures. The low crustal velocities are consistent 
with a hot or partially molten lower crust. 

7. Simple one-dimensional heat conduction calculations 
suggest that no anomalous mantle heat source is needed to 
explain the volcanic and tectonic activity provided that a 
heat flux of about 0.75 to 0.9 HFU (depending on the amount 
and distribution of crustal heat production) is available 

from the asthenosphere and that the chickness of the lithosphere 
does not increase drastically during crustal shortening. 
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on the second line for each event. 
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Figure Captions 

Fig. 1. Map of Tibet showing earthquake (dots) - station 

(triangles) combinations used for surface wave dispersion 
and the phase observation. Events I and II generated 
large long-period Rayleigh waves. Note the difference 
in the paths to stations from these two events. Event I 
was used for the Pl phase observation also. Major 
strike-slip and thrust faults (see legend) were also 
plotted following Molnar and Tapponnier [19751. 

Fig. 2. (a) Seismograuns of clear long period Rayleigh waves 

and PL-phase. These are the only Rayleigh wave observations 
found with a large signal to noise ratio at low frequencies, 
(b) Representative seismograms observed at stations south 
of the Himalaya for earthquakes in Tibet. Note that no 
clear Rayleigh wave signals with period beyond about 30 sec 
were apparent. Love waves vary from pulse-like signals 
to dispersed wave trains up to about 60 sec period. 

Fig. 3. (a) Group velocity dispersion curves of Rayleigh 

waves. No clear signal was observed between 30 to 60 
sec periods. A 3% maximum error in group velocities 
would be about 0.1 km/s. (b) Group velocity dispersion 
curves of Love waves. 
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Fig. 4. Phase velocity dispersion curves for Rayleigh waves. 
Observed values for I-NDI path (crosses) and II-SHL 
path (dots) were plotted together with three theoretical 
curves calculated for models F4/ S7, and S4. Model SE 
generates nearly identical phase velocity dispersion 
curve as that of the model S7 at the period range 
shown. Error bar shown is approximately +2.4% maximum 
error in phase velocity. See Figs. 5, 6, and 9 for 
model parameters. 

Fig. 5. Theoretical Rayleigh wave group velocity dispersion 
curves (top) for three models (bottom) with a crustal 
thickness of 55 km for different uppermost mantle shear 
wave velocities (Bj^) . The two long-period observations 
are shown as solid triangles for comparison. Also shown 
is the theoretical Love wave dispersion curve for model F4. 

Fig. 6. Theoretical Rayleigh wave group velocity dispersion 
curves (top) for 3 models (bottom) with a crustal 
thickness of 70 km for different 8^. Layout sams as 
in Fig. 5. 

Fig. 7. Theoretical Rayleigh wave group velocity dispersion 
curves (top) for 3 models (bottom) with a crustal 
thickness of 85 km for different Bj^. 
refered to. Layout same as in Fig. 5. 


Model E4 is not 
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Fig. 8. (a) Theoretical Rayleigh wave group velocity 

dispersion curves (top) for 3 models (bottom) with 
■ 4.4 Icm/s and a crust of 53 1cm thick but with 
different velocity structures in the crust. The average 
crustal shear velocity ( 7 ) ranv;;es from 3.49 km/s to 
3.44 Icir./a. The model FF3 which fits the observations is 
the same as model F4 (Fig. 5) . The upper and lower bound 
of our observations are shown in reversed solid triangles 
and solid triangles respectively. (b) Theoretical Rayleigh 
wave group velocity dispersion curves (top) for 4 models 
(bottom) with a fixed crustal thickness (55 km), (4.4 
km/s) and ^ (3.44 km/s). Note that a uniform crust 
between the Moho and the sediments with this low does 
not fit the observations. 

Fig. 9. (a) Theoretical Rayleigh wave group velocity dispersion 

curves (top) for 3 models (bottom) with fixed crustal 
thickness (70 Icm) , and (4.7 km/s). ranges from 
3.53 to 3.45 )cm/s. For all cases the velocities are too 
high. Layout same as Fig. 8. (b) Theoretical Rayleigh 

wave group velocity dispersion curves (top) for 3 models 
(bottom) . Model SS3 in (a) is plotted for comparison. 

The addition of a relatively high velocity layer 
(S a 3.9 km/s) makes model SS5 (= model S7, Fig. 6) 
fit the data at both long period and short period ends. 
Model SE is described in (c) . (c) Theoretical Rayleigh 

wave group velocity dispersion curves for models SE and 
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SSS. The slight increase in in SE has the same effect 
as the introduction of a layer at the base of the crust 
with B ■ 3.9 km/s in SSS (see (b) ) . Also note that 
model S7 (SSS) and SE have nearly identical phase velocity 
dispersion curves (Fig. 4). 

Fig. 10. (a) Two synthetic seismograms for different Q 

structures together with low-pass filtered vertical 
component Rayleigh wave of the II-SHL path observation. 
Both synthetics employ model S7 (Fig. 6) with source 
depth S km and an arbitrary seismic moment. All the 
seismograms have been normalized to have the same 
maximum peak-to-peak 2 unplitude» emphasizing their 
similarity in wave forms, and for each the portion 
earlier than 590 sec has been truncated. (b) Three 
low-pass filtered synthetic seismograms of the vertical 
component of Rayleigh waves for model S7 with different 
source depths of 5, 10 and 30 km. Amplitude scales havo 
been normalized as in (a) . (c) Low-pass filtered vertical 

component of the Rayleigh wave signal from the I-NDI 
path together with three synthetic seismograms for 
different velocity structures: models F7 , S4, S7 
(Figs. 5 and 6). Amplitude scales have been normalized 
as in (a). Note the particularly poor fit of model F7 . 
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Fig. II. (a) Vartlcal component of the long-period phase 
for the I-NDI path (solid trace) and one synthetic 
oeismogram (dashed trace) for a 70 )cm thiclc crust 
(single layer over a half-space) model, with an arbitrary 
seismic moment. The large amplitude P wave at the 
beginning of the observation results from deep mantle 
reflections and refractions that were not synthesized. 

(b) Two synthesized vertical components of Pl phases for 
layer-over-half space models with different crustal 
thic)cnesses of 70 and 35 ]<m. P and S wave velocities 
for crust and mantle are shown in the lower-right corner. 
Note the similarity between the two waveforms. (c) 
Synthesized vertical component of a P^^ phase for model 
SEMI (dashed trace) together with the observed signal. 

The inclusion of a crustal reflector alters the long- 
period waveform considerably from the pulse li)ce 
character in (b) . 

Fig. 12. Plots of Pn travel times normalized for a focal 

depth of 33 )cm (top) and travel time residuals (bottom) 
vs. distance to LHA for earthqua)ces north of the 
Indus-Tsang Po suture. The distance range is about 
3®-17®. The velocity (8.15+0.04 )cm/s) and travel 
time residuals were obtained by a least-square fitting 
of the 61 observations assuming a Moho depth of 70 1cm. 
Most of the residuals fall within +1 sec of the mean. 
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The dashed line shows the slope of a curve to which 
the data would be parallel if the velc;ity were 
7.9 km/s. Travel time residuals can be viewed as 
reduced travel times in conventional refraction analysis. 

Fig. 13. (a) Plots of P^ travel time residuals without depth 

normalization vs. distance to LHA for earthquakes north 
of the Indus-Tsang Po suture. Distance range is B^-IO**. 
(b) The same dataset as in (a) but with a normalization 
for a focal depth of 33 km. Note the reduction of 
residuals comp^u:ed with (a). Layout as in Fig. 12. 

Fig. 14. Plots of Pf^ travel time residuals normalized for a 
focal depth of 33 km vs. distance to LHA for earthquakes 
south of the Indus-Tsang Po suture. Distances range from 
3 « to 10®. Layout same as in Fig. 12. 

Fig. 15. Plots of Sj^ travel time residuals normalized for 

focal depths of 33 km vs. distance to LHA for earthquakes 
north of the Indus-Tsang Po suture. Distances range 
3® to 10®. Layout as in Fig. 12. 

Fig. 16. Map of Tibet showing the location of LHA and the 
epicenters used for the refraction profiles. Tect lic 
features are the same as in Fig. 1. 
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Fig. 17. Map of Tibet and the Himalaya ahowing the locations 
of 17 earthquakes (circles) used in teleseiamic travel 
time measurements. Roman numerals inside the circles 
are the S-P travel time residuals with respect to the 
J**B tables. All the values are rounded to the smaller 
nearest integer and are given at 1 eec intervals. 

Source depths have been normallzod to 5 km for earthquakes 
in Tibet and to 33 km for those in the Himalaya. Tectonic 
features are the same as in Fig. 1. ' 

Fig. 18. Examples of seismograms used in the teleseismic 
travel time measurements for earthquakes in Tibet and 
the Himalaya. Arrows indicate the onset of S wave 
arrivals. Numbers in the parentheses are S wave travel 
time residuals in sec with respect to the J-B table for 
the source depths given by the ISC. 

Fig. 19. (a) Plots of teleseismic P and S wave travel time 

residuals with respect to the J~B tables vs. distance 
for earthquakes in Tibet (b) S-P travel time residuals 
(circles) . The crosses are the results after station 
corrections have been added. Dashed line shows the 
zero-residual level for a focal depth normalized to 
5 km. The averaged S-P delays are given in Table III 
with estimated standard errors. 
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Fig. 20. (a) Plots of taleseismic P and S wav* travel time 

residuals with respect to the J-B tables vs. distance 
for earthquakes in the Himalaya. (b) S-P travel time 
residuals. Symbols are the same as in Fig. 19. 

Fig. 21. Average S-P travel time residuals for earthquakes 
in Tibet (left) and the Himalaya (right). Hatched and 
open columns are average residuals without and with the 
station corrections, respectively. The residuals are 
given for three different focal depths: 5 )cm (top), 

33 km (middle) and the depth reported by the ISC (bottom) . 

Fig. 22. Solidi of selected possible crustal and upper mantle 
materials of different water content from data obtained 
from experimental petrology (Lambert and Wyllie, 1970, 
Merrill et al., 1970, Merrill and Wyllie, 19751. Dashed 
lines indicate areas of large uncertainty cited in the 
original literature. 
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Abitroct. Th# upporwot aontlt F wbv# 
voloeltltB btnoBth Turkoy tod Iron voro ••tlaotod 
by Bpplylnc convontlonal trovoi tlM-dlftABCo 
rolBtion Mthod to orrivol tlMt of woll locBttd 
•orthquBkoo roeordtd Bt b fow itBtlonB. Tbo 
BvorocB upponoit BBntlo ? vbv# voloelty undor 
Turkoy it ottiiBtod froa two ttBtloot of tho 
World Wldt StBBdBrdlxod SoliBosrBph Nttvork 
(VWSSR), IttBBbul tnd T.^brif. Tho dBtB bt* 
eoniittont with » erutt of unlfoni, but poorly 
dotoralntd, thleknttt BBd bb upporaott aantlt P 
wBvo voloelty of 7.T3 ♦ 0.08 ka/t. Thlt voloelty 
It vory tlallBr to that for tho Ao«obb Sob BBd 
tUMOttt thBt Itt Btrueturo could bo elotoly ro- 
iBtodtothBt bOBOBth Turkoy. For Iron, tho 
roBultt cBlculBtod froa trovol tlaot to throo 
WV88R ttBtloBt, Mothod, Shlroi BBd TBbrlz. cbb 
bo oxplBlBod by t erutt diralBS toword tho 
touth-touthoBtt Bt Bbout 1* with bb upporaott 
BBBtlo P WBVO voloelty of 6.0 0.1 ka/t. If 

tho eruttBl thleknott woro 3k la 1 b tho Borth It 
would roBch Bbout kp ka 1b tho touth. Batod oo 
thoBo upporaott aoBtlo voloeltlot, tho toaporo- 
turo Bt Hoho benoBth Turkoy It probably eloto t>) 
tho aoltlBg toaporaturo of porldotlto but that 
boBOBth IrBB It probably lowor. 


iBtroductlon 

Tho BBBtlo voloelty and eruttal thlckoott aro 
iaportBBt eoattralntt 1 b UBdorttaBdlBC tho 
toetonlc ovolutloB of bb aroa. Thoro aro no 
rollablo rulot for prodletlof thoao Inforaatlona 
fron turfaeo olovatiOB or othar goophytlcal 
data. Tho Batin and Ran«o provlneo and Tlbot 
art good ox^plot. Both aroat havo high turfaca 
olovatlon, active faulting and roeont volcanlea 
(o.g.. 1976; Prlottloy 

and Bruno . 1976] . With only turfaeo vavo 
dltportiofi, tho Batin and Raago provlneo vat 
oneo thought to havo a thick erutt l Ewlna and 
Pratt . 1959]. Lator, dotallod rofractlon and 
additional turfaeo vavo dltportion atudloa of 
cruatal and uppor aantlo atrueturo In thla aroa 
[o.g., PraAtta. 1970i Prlottloy and Bruno. 1978] 
thowod that tho erutt la thin (^30 ka) and that 
tbo upporaott aantlo P vavo voloelty It low 
(7.8 ka/t). Tlbot, on tbo otbor hand, hat a 
thick en^ft (TO ka) and tho upporaott aantlo P 
vavo voloelty (8.1 * 0.1 ka/t) It alaoat tho 
tamo at that of ttablo ahlold aroat ( Chon and 
Molnar , I960]. 

Turkoy and Iran llo north of or vlthln tho 
ccntlnontal collltlon zono botvoon Africa and 
Guratla. Tho turfaeo foaturot in thoto tvo aroat 
aro vory tiailar to tho Batin and Range provlneo 
and Tibetan region. They all havo active 
faulting, roeont volcanlea and high turfaeo 
olovatlona. In tbit ttudy, tho uppomott mantle 
P \ VO (Pq) voloeltlot under Turkoy and Iran aro 
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Invoatlgatod and coaparod with tho tvo provlncot 
■ontlonod above. Bated on thoao ottlaatod Pn 
voloeltlot, VO place conttralntt on tho upper- 
■OBt aantlo toaporaturot beneath Turkoy and Iran. 

Data Analyalt 

Vo oxaalnod all tbo oarthquakot that occurred 
In Turkoy and Iran during January 196 k to 
January 1976 and for which travel tlaot woro 
reported at latanbul (1ST), Tabriz (TAB), 

Mothod (NBH), and Sblrat (8HX). Tho location* 
of tho oartbquakea uaod In thla ttudy aro thovn 
In Fig. 1. Tho origin tlaot, focal doptht, 
opleontral dlttaneot and arrival tlaot woro 
taken Trcm tho Bulletin of tho International 
Soltaologleal Centro. Tho oarthquakot wore 
aortod baaed on tbo following criteria: (1) wo 
uaod only oventt for which aoro than 25 
aoaauroaoDta of P vavo arrival tlaot voro utod 
In tho dotoralnatlon of tho earthquake bypo- 
eontor and origin tlao; (2) oventt with reported 
Btandard dovlatlona of tho origin tlaot greater 
than tvo aoeonda voro oxeludodi (3) only oventt 
with focal doptht lota than 35 ka wore Included. 
Initially wo eonaldorod oarthquakot with 
reported foeal doptht at largo at 60 ka, but 
data teattored very aueh for thoao deeper than 
35 ka. So, only oventt with foeal depth lota 
than 39 ka ware Ineludodt (k) only *vontB wltl 
opleontral dlttaneot botvoon k^ and 11. 9° woro 
Included. Within tbo opleontral dlataneo range 
of approxlaatoly 3° to 12°, the flrtt arrival 
on tho tolaaograa for a given event it the P^ 
phaao. Tho exact range depondo on tbo cruatal 
thleknott and tbo P vavo voloeltlot In the erutt 
and aantlo. Epieontral dlttaneot botvoon 3° 
and 12° wore axaalBed for all atatlona, and 
alaott all the travel tlao rotldualt with 
rotpoet to tho aoan aro negative between 3° and 
k°. Alto, oventt oeeurred In Iran with 
opleontral diataacot at largo at lk° reported 
at Tabriz wore Included and the travel tlao 
rotldualt voro all negative for oventt with 
dlttaneot greater than 11.9°. Tho oventt with 
dlttaneot lota than k° or greater than 11 . 5 ° 
woro excluded for all tbo four atations to 
aalntaln a tingle dlataneo range for th* entire 
data tot. 

After applying that* eritoria, there wore 
Btlll toB* oventt with largo rotldualt (greater 
than two toeondt) at on* or nor* atatlont. Upon 
chocking th* Boitaograaa w* found that for mott 
of that* eaaoB tho arrival* wore oaorgont. The 
algnal to nola* ratioa are vory low and th* 
arrival tlaot could not b* dotorslnod precltoly 
for aany of thea. Thoao ovonta wore excluded 
froa thla ttudy. Sob* of the reported arrival 
tlaot war* rovltod after chocking tho toltmo- 
graaa. Fig. 2 abeva an oxaapl*. 

For a particular oarthquako 1, let T^ be the 
travel tlao calculated froa th* obtorvou arrival 
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ri«. 1. Nftp ihowlnc tpiotntitrt of oorthquokot uaod for •■tlattinf tht upporaott winUt P wavt 
voloeltlts bmootb Tufkoy and Iran. Cirolaa, orosaoi. pluaoi. and aquaraa show tarthquakaa 
vitb arrival tisoa raportod at latanbul, Shlrat, Tabrix, and Ntahad, rtapactivtly . 


TABLE 1. Suaaary of P„ Valoclty, Inttretpt TIm and Crustal Thlcknoss Undsr 
Turksy nnd Iran Estlaatad Proa Arrival Tlaas to Dlfftrsnt Stations. 



Recorded 

Muabers of 

Before dspt 

h-corrsctlQfl 

After dsbth-corrsctlon 

Region 

station 

earthquakes P„(ka/s) Splssc) F„(ka/s) ap(csc) 

iBcludsd 

crustal 

thickness 

Turkey 

nf 

TAB 

ff 

23 

7.60X0.06 

7.01X0.08 

3.3X5. 3 

5.8x'>.? 

7.70x0. Ok 1. 7X3.1 
7.76X0.06 3. 9X3. 9 

25x32 km 
36x3^ km 

Iran 

MSB 

SMI 

TAB 

39 

50 



9.02x0.05 

9.19x0.15 

7,79x0.06 

7.kt0.£ 

8.7X13.1 

8.02X0.07 6. 5X6. 6 

e.lktO.lO 7.1x6.k 

«7xke ka 
k9x55 km 
3‘*1‘»9 km 


tlas and tbs rsportsd origin tla« yotalnsd by 
usinf aany stations in a global seals. Ws 
dsflns tns travsl tlas rssldual by 

Ci • Ti . (7^ ♦ ap). 

Tbs uppsraost aantls P wavs vsloelty (V-) and 
Intsrespt tlas (sp) wars ealculatsd by Isast- 
squara fitting or travsl tiass (T|) vsrsus 
dlstauca (bi) and ainlalsing tbs sisi of ths 
squarss of tbs travsl tlas residuals (lc<^) for 
aartbquakas rsportsd at tbs particular station. 

For eoaparlson, both dspth-eorrsctsd travsl 
tlnss and travsl tiass vltbout dspth corrsctlona 
wars ealculatsd. For tbs dsptb-eorrscted case 
all of tbs origin tiass and tbsrsfors tbs travsl 
tiass, wars nonallssd to a focal depth of 33 ka, 
which Is tbs asan rsportsd focal depth. To do 
this, vs assuBsd a crustal P wavs vsloelty of 
6.1 ka/s and aantls P wavs velocity of 9.0 ka/s. 
Different ass\aptlons of tbs velocity structure 
will not give large differences in tbs results. 

For sxaapls, asstalng Pn of 7.7 ka/s or 3.k 
ka/s for sn earthquake occurring at tbs sartn's 
surface aakos only O.k seconds difference In 
ths calculated dsptb-corrsctsd origin tlae. 

Results 

Ws used 37 earthquakes rsportsd at 1ST and 23 
earthquakes rsportsd at TAB to calculate ths 
average uppensost aantls P wavs velocity bsnsatb 
Turksy. For ths eases without a dspth-eorrsetioa. 


tbs Pq vslucltlsa are estlaatad to be 7.60 ♦ 0.06 
ks\/' and 7.81 s O.06 ka/s froa data reported by 
IPT and TAB, respectively (Table 1). For the 
ds^ ch-eorrec ted cases, the Pn velocities ere 
s .iaatod to be 7.70 ♦ 0.0k ka/s and 7.76 ♦ 0.06 
ka's respectively. Froa Fig. 3 and the earculsted 
standard deviations, we can see that the scatter 
of the data Is reduced by the depth>correctlon 
as one aay expect theoretically. Figs. 3d and 
3e show that If the velocity were forced to be 
6.0 ka/s, the data should define a line with e 
slrpe parallel to ths dashed lines. Froa these 
figures, we can see that a high ?„ velocity will 
not fit the oata wall unless there is soae 
distance-dependent systsaatlc error In dater- 
Binlng tbs spleantars, focal depths and origin 
times . 

The large standard deviations In ths Intercept 
tines (Table 1), sake It difficult to estimate 
ths crustal thickness froa these values. Since 



Fig. 2. Selsaograa shows the first arrival 
of an event reporte<' at Istanbul, Aug. 17, 

1971 (0l*hr29Bin33.usec. 37,07®S. 36.77®E'. 

Dashed arrow indicates the reported reading 
(0khr3iBln2i.fl**C) , and solid arrow indicafes 
ths first arrival (Okbr3xaini9seo ) that we used. 




ChM tt al.: f lava Valocitiaa lanaath Itirtay and Iran 


79 


^200^r 


I 


l«O.Oh 


1 20.0 H 


80.0 K 


S 40.0 1 



U--i- —I- I l_ — I 1— 

T5 £5 a5 ioo 

DISTANCE (DEG) 


12.0 





Pig. 3. Plots of travel tloaa and travel tine 
residuals vs. distance, (a) and (b) show 37 
travel tines and residuals without correction 
for different focal depths rscorded at X8T, and 
(c) and (d) show these tines corrected for 
different focal depths, (e) shows travel tine 
residuals corrected for different depths of 
focus for 23 earthquakes In Turkey reported at 
Tabrls. Sash lines In (d) and (e) show the 
slops of a line to which the data points would 
be parallel If the Pq velocity were 8.0 kn/s. 


no significant differences appear in the Pn 
velocities and Intercept tines estlnated fron 
arrival tines reported by 1ST and TAB, the 
aasuqttion of a flat layered crust for Turkey 

seens reasonable. 

The velocity beneath Iran, estlnated fron 
three stations, MSB, SHI, and TAB are 

8.02 s 0.09 kn/s, 6.19 *. 0.15 kn/s and 
7.79 ± 0.06 ha/s, respeFtively, for the eases 
without m t..jpth correction, and 6.02 s 0.07 
kn/s, 6.1b ♦ 0.10 kn/s and 7.85 *, 0.05 kn/s. 
respactlvely, for the depth-corrected cases 
(Table 1). These data could be explained by a 
sinple nodal of a crust dipping towards the 
south-southeast at about 1° with the thicker 
part near Shiras and the thinner part near 
Tabriz. The average uppemost nantle P wave 
velocity would bo about 8.0 ♦ 0.1 kn/s. Again, 
the uncertainties of the Intercept tines (Table 
1) are very large, but this could be partly due 
to the dipping Moho. 

Considering the scatter of the data (Fig. b) 
allowing large uncertainties for those velocity 
estlnatos, it is possible that the crustal 
thickness Is roughly unlfons (bo-50 kn) beneath 
Iran with nantle velocity around 8.0 kn/s. 
niase ostlnates of the crustal thickness and Pn 
velocity are average values and nore conpllcated 
interpretations are also possible. 
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Discussions and Inpllcations 

The Pjj velocity under Turkey ( 7.73 ♦ o. 08 
kn/s) Is conparable to that of the BasTn and 
Bonge province (7.8 ^ 0.1 kn/s) fa.g« Prodehl . 
1970 ] and is uigaifi?antly lover than that under 
Tibet (6.1 ♦ 0.1 kn/s) ( Chen and Molni^ . I960]. 
The estlnated P^ velocity under Iran (8.0 * 0.1 
kn/s) is higher than that under the Basin and 
Range province and lover than that under Tibet. 
The Pn velocity under Turkey and Iran differ 
significantly fron each other although the two 
regions are adjacent to each other and both lie 


•0 
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Fig. U. Plots of traval time residuals 
corrected for different focal depths vs. 
distance for earthquakes in Iran reported at 
Meshed (a). Shiraz (h), and Tabriz (c), 
respectively. Layout as In Figure 3. 


north of or within the continental collision 
zone extended between Afro>Arabia and Eurasia. 

The Pq velocity beneath Turkey Is very slailar 
to that beneath the Aegean Sea (7. 6-7. 8 km/s) 
froB deep selsale soiadlng [ Makrls . 1978]. This 
suggests that the aantle structure beneath the 
Aegean Sea could be closely related to that 
beneath Turkey [ McKenzie . 1978]. 

With these velocltiea, we can place rough 
constraints on the tea^eratures at the Moho 
beneath Iran and Turk^. Assuming that the 
temperature and pressure dependencies of the 
P wave velocity In the \q>perBost mantle are 
similar to those of the, olivine, then a pressure 
decrease of 10 kb or te^rature increase of 
250°C will reduce the Pn velocity by 0.1 km/s 
[e.g., Birch . 1969]. 

The temperatiire In the uppezvost mantle 
beneath the shield areas estimated from heat 
flow data f Sclater et al. . 1979] If about UOO® 

♦ 100°C. Taking this temperature and average 
7hleld Pq velocity of 8.1 km/s and average shield 


crustal thickness of 35 km, with Pq velocity of 
8.0 e 0.1 km/s, the uppermost mantle temperature 
beneath Iran would be 650" e 350"C if the 
crustal thickness were 35 km, or 800 "r 3$0"C 
if the crustal thickness were 55 km. ~Wlth Pq 
velocity of 7.73 ± 0.08 km/s and crustal thick- 
ness of 25-35 km, however, the temperature In the 
uppermost mantle beneath Turkey would be 
calculated to be IbOO" ♦ 300"C. This tempera- 
ture would be high enouih to melt perldotlte at 
uppenost mantle pressures. Hear the melting 
temperature, the asaiaed linear relationship 
between the Pn velocity and tnperature at 
constant pressure le likely to be invalid 
[ Qoetze. 1977]. Indeed, the actual temperature 
in the uppenost aantle beneath Turkey is 
surely lower than lkOO"C, but probably is above 
the perldotlte solidus. This estimate is 
consistent with the recent basaltic volcanlsm 
in Turkey and implies that the basaltic liquid 
originates la the uppermost mantle. 
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ABSTRACT 

Haui, A.T. and Tokadi, M.N., 1070. The evolution of thermal atructurea beneath a aub* 

duction aonu. Tectonophyaica, 60; 43—60. 

Thermal modela of the aubducting lithoaphere and the aurrounding mantle are inveati- 
gated aimultaneoualy. In the model, alab aubductlon ia treated u a kinematic proceaa. 
Subduction angle and aubductlon velocity are preMribed. The mantle ia permitted to 
reapond dynamically to aubductlon. Mantle flow ia driven by both mechanical ahearing 
and thermal buoyancy. Finite d'ff.*’ence numerical techniquea are utilized to model the 
subduction process. Results indicar sat the descending lithoaphere and the induced flow 
in the wedge area dictate the the *. al structure beneath plate convergent regions. The 
temperature of a slab Is found to be insensitive to the mantle condition prior to slab sub- 
duction. Viscous shear heating effects within the dynamic mantle are found to be impor- 
tant in maintaining a high temperature region immediately beneath the inter-arc basins. 


INTRODUCTION 

The subduction regions are among the most complex tectonic provinces 
on Earth. Many processes taking place in these regions have not yet been 
understood satisfactorily (Uyeda, 1977). The thermal regime of the region is 
undoubtedly a piece of important information which enable us to better 
understand many geophysical and geological processes taking place in the 
converging plate boundaries. The thermal structure of a downgoing slab has 
been atudied by a number of investipstors (McKenzie, 1969; Turcotte and 
Oxburgh, 1969; Minear and Toksdz, 1970a, b; Hasebe et al., 1970; Toksoz 
et al., 1971, 1973; Turcotte and Schubert, 1971, 1973; Griggs, 1972; 
Schubert et al., 1975). These models were calculated without taking into 
account the effects of the dynamic mantle. Using these slab models as bound- 
ary conditions, the thermid structure of the surrounding mantle has also 
been investigated (Andrews and Sleep, 1974; Jones, 1977;Toksbz and Bird, 
1977; Toksdz and Hsui, 1978). These two families of models generally 
ignored the continuous interaction between the slab and the mantle. In this 
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pap«r w« conitruct • mort raalUtic model beneath the lubduction region by 
calculating the temperature within a liab and that of the (unrounding dy- 
namic mantle limultaneoualy. 

Another problem we investigate is why the temperature above the slab is 
relatively high. High temperatun in the wedge area is inferred by the high 
surface heat flow meuurements (for a good summary, see Watanabe et al., 
1977) and the high attenuation of seismic waves (both P and S) in these 
regions (Oliver and (sacks, 1967; Molnar and Oliver, 1969; Utsu, 1971; 
Barazangi and (sacks, 1971; Barazangi et al., 1975). Beneath the converging 
plate boundaries, the mantle is cooled by the subducting oceanic lithosphere, 
(f the subduction rate is lOcm/yr and the plate thickness is 100 km. a volume 
of about 10* km’ of cold lithospheric material will be injected into the mantle 
per million years per km width of the subduction zone. Taking an average 
temperature difference of about 500* C between the lithosphere and the 
mantle beneath, it follows that the same volume of the mantle will be 
replaced by the 500*dcpee cooler lithospheric material in a million years. 
Within the same period of time, radiogenic heat production (assuming that 
heat production is at equilibrium with the surface heat fluxes and the heat 
sources are uniformly distributed within the mantle) is only able to increase 
the temperature by about 1*C. Then, what is the mechanism that keeps the 
wedge region from cooling? This paper looks at this problem in some detail. 

MODEL DESCRIPTION 

Ideally, it would be desirable to treat the subduction and mantle-convec- 
tion problem as a closed, continuous and dynamic system with appropriate 
mantle properties. However, this presents extreme difficulties, both analyti- 
cally and numerically. In this paper, the thermal structure of the subduction 
region is studied by modeling a dab subducting kinematically into a deform- 
able mantle. In other words, the subduction velocity and the angle of sub- 
duction are prescribed. Instead of posing the problem as a continuously 
evolving system, we calculate the thermal regime at discrete time intervals. 
The basic approach is very similar in philosophy to that described by Mineur 
and Toksdz (1970a, b) such that we are in essence studying the heating of a 
discretely descending slab by the dynamic mantle. The slab is moved down- 
ward by one grid space at a time. Calculations are then carried out to deter- 
mine how the slab is heated up by the surrounding mantle within the time 
interval necessary for the slab to reach that depth. The convection and the 
thermal regime of the mantle are also determined simultaneously. Our model 
is an extension of most of the published slab models. 

The model consists of a 100 km thick conducting layer (lithosphere) lying 
on top of a convecting 700 km thick upper mantle. Calculations for our 
model are separated into two parts. One part is to calculate the thermal 
structure within the lithosphere and the subducting slab, while the other is 
to compute the thermal regime within the dynamic mantle that lies beneath 


th« conducting Uthotphcrt and nirrounda the dc'ctnding slab. The two cal- 
culationa interact with each odier through their common boundaiiei. When 
the dab subducts, its tMnperature is translated to a new position, the mantle 
is then allowed to respond to this disturbance within the time interval neces- 
sary for the slab to travel to its new position. The effects of slab shearing, 
thermal buoyancy, and adiabatic compression are incorporated in the calcu- 
lations. Computaitions of the slab temperatures and the mantle temperatures 
are self-consistent at each time step through the matching of boundary 
conditions. An explicit, forward time march, finite difference numerical 
scheme is used to simulate tiie thermal consequences of the subduction 
process. Details of the mathematical formulation, the numerical approxima- 
tions, and the computational procedures are very similar to those given in 
Toksbz and Hsui (1978). 

Boundary conditiom. The boundary conditions are difficult to specify since 
neither the mantle flow structure nor the nature of the bottom of the mantle 
convection layer is known with great confidence. However, if the vertical 
boundaries are sufficiently far away from the subduction region, the effects 
of the vertical boundary conditions will be small on the thermal structure of 
the slab and its immediate environment. In our models, the horizontal 
dimension spans about 2800 km. The location of subduction at the surface is J 

chosen at least 800 km away from either end, which is more than the vertical 
length scale of the convecting mantle. After the slab has been translated to 
its new position, the thermomechanical structure of the mantle is calculated 
assuming the vertical boundaries are shear streu free, mass exchange free, 
and thermally insulated. At tiie top of the dynamic mantle, a prescribed 
velocity equivalent to the slab velocity is imposed on the oceanic side of the 
slab while a no-slip condition is imposed on the continental side. The outward 
heat flux from the mantle at this boundary is determined from the thermal 
structure of the overlying lithosphere. At the bottom, a prescribed inward 
heat flux is impMed and a no-slip mechanical condition is used. We also cal- 
culated models with a shear-etress-free bottom boundary condition. The ther- 
mal structures are not very different from those of a no-slip bottom boundary 
condition. This is in agreement with the conclusions of Hsui (1978) that 
thermal convection solutions are not sensitive to the mechanical boundary 
conditions at the bottom. The most difflcult boundary condition to choose 
is probably that at the tip of the slab. In this paper, we choose a no flow 
condition at this boundary, so that it is consistent with the no mass flux con- 
dition imposed on the two vertical boundaries. Mass flow induced by the tip 
strongly depends on the shape of the tip. Thus the choice of this boundary 
condition may introduce different effects locally. Richter (1977) studied the 
effects of tip push and found that it could be siftnificant. Howevet the 
effects of variable mantle viscosity and different boundary geometries are 
not investigated. Furthermore, when one introduces a mass flux at the tip of 
a slab (mass inpi.t) one also has to specify the distriburion of outward 
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maas flux at other boundaries. Different distributions of outward mass flux 
will affect the solutions in their own way. Therefore, this alternative will 
also introduce its own uncertainties to the problem. It should be stated that 
the “tip** problem goes away when the slab penetrates the bottom. 

initial eondition$. Since the problem is cast as a bransient problem, initial 
condititms are required to complete the mathematical formulation. Again, 
neither the mechanical nor the thermal conditions within the mantle prior to 
slab subduction arc known. Therefore, we choose to treat the initial condi- 
tion u a parameter and study the sensitivity of the tftermal regime of a slab 
and its immediate environment to the initial conditions chosen. Even though 
very little data exists which can be used to constrain the initial conditions, 
we auume that the mantle was probably at equilibrium with the heat inputs 
to the system prior to slab subduction. Therefore, a steady state solution for 
the dynamic mantle in the absence of a slab is obtain^ and used u the 
initial condition for the mantle. Generally, any solution of thermal convec- 
tion will yield upwelling as well as downgoing currents. In order to study the 
local effects upon a subducting slab, subduction has been initiated at three 
different locations: above the upwelling current, above the downgoing cur- 
rent, and above the center of a convection cell. With these three modes of 
subduction, the regional effects of different initid conditions on slab tem- 
peratures can be analyzed. 

Heat tourcee. The mantle not only deforms in response to the mechanical 
shearing of a descending slab, but motions can also occur within the mantle 
due to thermal buoyancy generated by internal heat release. Therefore, the 
specification of heat source strength within the system is necessary. How- 
ever, the distributions of U, Th, and K heat sources within the Earth are not 
known in detail. In this paper, the strength of the heat production rate is 
estimated by assuming that it is at equilibrium with the surface heat flux. 
The average heat flow values from the ocean floor vary from 50—80 erg/cm ‘-s 
(Sclater et al., 1979). We further assume that the lithosphere and the lower 
mantle >«ach contribute about 25% of this heat budget while the rest is 
generated internally within the upper mantle, (i.e.. The heat production rate 
within the lithosphere is about 1.5 ■ 10** erg/em’-s. The heat production rate 
within the upper mantle is 4.5 * 10*'’ erg/em’-s and the heat flux from the 
lower mantle is about 13.5 erg/em’-s.) The choice of heat source distribution 
does not appear to affect the results very much since the time duration is 
short for the radiogenic isotopes to generate significant heating. For the 
same reason, the heat production rate is taken to be constant in time in our 
calculation. 

Viicoiity. The relationship between the rate of deformation and the shear 
stress within the mantle remains an obscure subject mamly because of the 
absence of relevant laboratory data at the high pressures and low strain rates 
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appropriate to the mantle conditions. However, numerical experiment! (Tur* 
cotte at al., 1973; Pannantier et al., 1970) indicate Uiat the thermal itruc- 
ture of a convective lyitem is very leniitive to the choice of a constant or a 
strongly temperature and pressure dependent viscosity, while relatively 
insensitive to the assumption of a newtonian or a non*newtonian behavior. 
Therefore, for mathematical simplicity, our model assumes a newtonian 
medium with a strongly temperature and pressure dependent viscosity. The 
viscosity v is given by: v" A exp(BT„/T), where T„ is the absolute melting 
temperature which is taken from the 0.1% hydrous pyrolite mantle model of 
Ringwood (1975). T is the local absolute temperature. A and B are two con- 
stants SO adjusted that the resultant viscosities are in agreement with those 
determined from the study of glacial uplift data (Cathles, 1975; Peltier and 
Andrews, 1976; Peltier, 1976). In this paper, A ■ 1.8 ■ 10'* cm'/s and S - 
9.0, are used. 

Pha$e boundariti. The upper mantle phase changes ( basal t*eclogite, pyro- 
lite-gamet, olivine-spinel, spinel-post spinel, etc.) could affect both the con- 
vection and thermaJ regime by density changes and heats of reaction. The 
effects of phase transitions on mantle dynamics have been studied by 
Schubert and Turcotte (1971) and Schubert et al. (1975). Their results indi- 
cate that the olivine-spinel phase transition it unlikely to inhibit mantle con- 
vection. Whether the spinel to post-spinel phase transition will inhibit mantle 
convection depends on the composition of the post-spinel phase. Since this 
transition it near the bottom of the convection zone, it 'viU have little effect 
on this problem. 

Besides the deformable mantle, phase transitions alto take place within 
the descending slab. Their effects on the thermal regime of a subducting 
lithosphere have been studied by Toktdz et al. (1971). They concluded that 
the contribution of phase transitions to the heating of the interior of a 
downgoing slab it small when compared with other heating mechanisms such 
at conductive and convective heat transport, adiabatic heating and shear- 
strain heating. Based on the latent heat associated with these phase transi- 
tions (Verhoogan, 1965; Akimoto and Fujisawa, 1968; and Ringwood, 
1970), the temperature increase due to phase changes within a slab it lest 
than 50* C locaUy. Therefore, in this study, we neglect the effects of phase 
transitions. 

THERMAL MODELS 

The numerical models presented here are obtained using a 57 x 15 (50 km 
grid spacing) x-r grid system for the convecting mantle and a 10 x 10 km 
grid system for calculating the temperature within the slab. A test cate using 
a 141 X 36 (20 km grid spacing) grid system for the mantle was carried out. 
Results are very similar to those obtained using the coarse grid. They differ 
by no more than 10 percent. Considering the uncertainties on many of the 
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Fig. 1. Th« initial Iharmal and mtchanical ttnictura of tha modal calculationi prataniad 
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*C ii givan by tha numbar aasociatad with aach iaotharm. Tha ihadad araa rapraaanta tha 
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1.24. Straam linaa plottad rapraaant lii'l/lii'initi ■ 1/3. 2/3. PoinUA, B. and C indicate 
tha locationa whara alab aubduction ia initiatad. 


parameters used in the models, the 50 km mantle grid system is adequate for 
calculating the gross thermal structure. 

Figure 1 shows the initial conditions we used for the models. The diagram 
on the top is a plot of isotherms (constant temperature lines), while the one 
at the bottom is a plot of stream lines. In the isotherm plot, the lithosphere 
(shaded layer) is superimposed on top of the convecting mantle. However, 
in the stream line plot, the lithosphere hu been removed. Points A, B, and 
C in the stream line plot indicate the locations where the initiation of sub* 
duction takes place for the three cases. In all the models, the slab subducts at 
an angle of 46“ counterclockwise with respect to the horizon and at a 
speed of 8 cm/yr. 

The following sequence of figures (Figs. 2—4) shows the thermal evolution 
history of a subducting slab and its surrounding mantle, for subduction initi- 
ating at point B. The shaded area represents the lithosphere and the subduct- 
ing slab. As in Fig. 1, the top diagrams in this sequence are plots of isotherms 
while the bottom ones are plots of stream lines. Figure 2 ^ows the thermal 
and flow structure at 5.0 m.y. after the initiation of subduction. The thermal 
fleld has not been disturbed greatly except near the descending slab where 
the gcotherms are being depressed dowrnward. However, the flow field has 
been altered signiflcantly. It is mainly because mechanical disturbances nor- 
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mally propagate much faater than thermal diaturbancei. Therefore, the 
whole region feels the existence of the descending slab and adjusts the flow 
pattern accordingly even though the slab has penetrated only a small 
distance into the mantle. As to the thermal structure, the low temperature 
geotherms are being carried downward into the mantle by the slab. In 
essence, this demonstrates the cooling of the mantle by the subducting 
oceanic lithosphere. As the lithosphere subducts, induced convection cur* 
rents are generated on both sides of the slab. Our models indicate that the 
induced flow above the slab (on the continental side and within the wedge 
area) has a higher intensity than that generated on the oceanic side of the 
slab. 

This phenomenon is not observed the flow is driven only by mechanical 
shearing (Richter, 1977; Tovish et al., 1978). Flow patterns driven by both 
mechanical shear and thermal buoyancy have been studied by Torrance et al. 
(1972). Even though the specific problems they studied are not related to 
the Earth’s mantle, the physical consequences they derived are directly appli- 
cable to our study. In their study, Torrance et al. (1972) concluded that 
buoyancy can significantly change or even dominate the flow driven by shear 
only, resulting in hi^er flow velocities. This is in agreement with our results. 
Due to tite combined effects of the large thermal gradient introduced by the 
descending slab, the strongly temperature dependent viscosity and a small 
flow cross section, the wedge flow bums out to be a strong current. Our cal* 



culAtion ihowt that th« flow voloeity in Uit wtdfe can rtach ai high at two 
or thraa timat tha ilab vaiocity. 

At 10 m.y. (Pif. 3), tha downfoinf ocaanic Uthoaphara hat raachad a 
dapth of about 550 km. Tha mantla it coolad on both lidat of the stab and 
tha inducad flow in tha wadga swmaint dominant. At 20 m.y. (Fig. 4). 
tha slab hat panatratad tha bottom of our modal. Evan though tha choica of 
a bottom at 800 km dapth it tomawhat arbitrary, compariton batwaan Fig. 
4 and Fig. 3 indicatai that tha tharmal ttmctura in tha immadiata environ< 
mant of tha tlab abova 300 km appaart not to be affactad tignificantly by 
tha pratanca of thit boundary. Sinca moat of tha gaophytical and gaologicai 
conatraintt lie abova 300 km dapth, tha attumption of a bottom at 800 km 
in our modalt will not affact vary much tha near nirfaca ratulta derived from 
thia ttudy. 

In tha ttraam line plota of Figt. 1—4, ttraam functiont are plotted at equal 
qiNKing for Mch diagram. Sinca tha maximum absolute ttraam functiont vary 
from one diagram to tha other at indicated in tha figure captions, tha dimin- 
ishing of ttraam lines does not imply the diminishing of the flow valocitiat. 
In fact, in all our calculations, flow velocities generally increase in time. 
Additionally, in our models, mantla viMosity varies from about 10'* poise 
imroadiataly beneath tha lithosphere to about 10’* poise near tha bottom. 
Tha viscoeity variation takas a vary steep gradient near the surface and a rela- 
tively shallow gradient away from the asthenosphare. However, the viscosity 
structure adjacent to the slab is very complicated. It can no longer be 
expressed in terms of depth only. Since viscosity is directly related to tem- 
perature, a complex viscosity structure around a slab is a direct reflection of 
the complex thermal structure around a descending lithosphere. 

Figures 5—7 raprMant t’ta mlargad view of the isotherm piots presented in 
the previous evolution sequence (Figs. 2—4). Only regions immediately sur- 
rounding the slab are included. Again, the shaded area indicates the litho- 
sphere and the descending slab. These diagrams provide better resolution of 
the thermal structure in the vicinity of a subducting lithosphere. In general, 
these diagrams demonstrate the cooling of the mantle u the low temperature 
isotherms are being carried downward by the descending slab. Since the man- 
tle on both sides of the slab is being cooM by the subducting cold lithosphere, 
the slab appears thicker thermaily except in the region between 100 and 200 
km depth. In this region, the cold slab is unnble to cool the mantle because a 
convective current is being induced in such a way that hot mantle material is 
continuously being pumped towards the vicinity of the slab at this depth. As 
a result, a bottle neck formation appears in the thermal structure of a slab 
s>.^bduction environment. Because of the coarse grid system used for numeri- 
cal computation, our models onnot resolve whether tiiis induced flow is 
able to maintain a temperature sufficiently high in the area immediately 
adjacent to the slab surface. This is important for the generation of magmas 
in island-arc volcanoes. Nevertheleu, our models demonstrate that the 
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Pig. 5. Enlargtd tamparaturt plot of Pig. 2 In tha vicinity of a daacanding ilab. 
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Pig. 6. Bnlargad tamparatura plot of Pig. 3 in tha vicinity of a daacanding slab. 
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Fig. 7. Enlargtd Umparatur* plot of Fig. 4 in tht vicinity of a dcscanding ilab. 


uthtnotphtro above the alab will not be golldifled by the deter n ling cold 
lithotphere, limply because the induced current will continuously feed hot 
mantle material into the area. 

The numerical results, showing that the temperature at the top of the 
wedge region remains relatively high, are consistent with the high attenua- 
tion of seismic waves and the hi^ surface heat flow of die back arc region 
(see flrst section for references). It should be pointed out that all the models 
presented thus far are constructed with die shear heating effect incorporated 
not only between the slab and the mantle, but within the mantle flow at 
well. If the shear heating effect within the mande flow is not included, the 
hifdi temperature region at the top of the wedge area no longer exists u 
shown in Fig. 8. Therefore, viscous shear heating within the mantle flow 
appears to be important in maintaining a hi^ temperature region beneath 
the back-ar? b?.«iiM. In fact, viscous ihenr heating can be u much as 2—3 
times higher than the radiogenic heating. Therefore, the requirement of shear 
heating to maintain a high temperature wedge does not appear too surpris- 
ing. 

In w.'der to see how much our slab temperature model differs from other 
models, a comparison is made (Fig. 9) between our model taken at 50 m.y. 
after the initiation of subduction and other pubUsl}ed models (McKenzie, 
1969; Toksdz et al., 1973; Turcotte and Schubert, 1973). The model by Mc- 
Kenzie (1969) represents the simplest slab model. It is analytic and considers 
only the effects of subduction and conduction. Therefore, it gives a lower 
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bound of the thormal rsfimc within ■ dt«c«nding iltb. Turcottc ind Schu* 
bwt (1973) havt alio dttarminad tha thannal itructure within a slab using 
analytic mtthods. The affect of phaae transition (olivine>ipinel) has been 
incorporated. However, because of the absence of internal heat generation 
and utilizing a lower temperature at the base of the lldiosphere, their models 
V>pear to be generally cooler than our models. Using finite-difference 
numerical techniques, Toksbz et al. (1973) have calculated the thermal struc- 
ture within a slab using a convective mantle geotherm and including the 
effects of adiabatic compression, internal heat release, shear heating and 
phase changes. This model appears warmer than our present model. Since 
our model has not considemd the effect of phase changes and because oli- 
vine-spinel phase change is an exothermic proceu, it can probably explain 
why our model is somewhat cooler. Meaningful comparison with the slab 
model of ToksOz et al. (1971) and Minear and Toksdz (1970a, b) cannot be 
made, however, because these models used conductive mantle geotherms 
that result in higher mantle temperatures. The present slab model is 
bracketed by other models u can be seen in Fig. 9, and it represents a good 
average temperature model for the subducting lithosphere. 

OTHER MODELS 

As pointed out previously, the condition of the mantle prior to slab sub- 
duction is subject to great uncertainty. Whether the thermal structure of the 
descending oceanic lithosphere and its environment are related to the initial 
mantle conditions is not immediately apparent. In order to investigate the 
relationship between the initial mantle conditions and the subsequent ther- 
mal evolution of a descending slab, we carried out models for which subdue- 
tion has been initiated at different locations. Figures 10 and 11 show the 
thermal structure of a slab and its environment at two different times. In this 
case, subduction has been initiated at a point where downgoing mantle flow 
is taking place (point A in Fig. 1). Comparing Figs. 10 and 11 to Figs. 5 and 
7, respectively, the thermal structures within the ilab and within the wedge 
region are found to be very similar. When subduction is initiated at point C 
(in Fig. 1) where upwelling mantle flow is taking place, results (Figs. 12 and 
13), again, are found to be similar to those where subduction is initiated 
either above a downgoing mantle flow or at the middle of a convection cell. 
From tiiese results, it is evident that the thermal structure beneath the sub- 
duction area is controlled by the descending lithosphere and the resulting 
induced flow in the wedge axm. The subduction process dominates whatever 
the mantle condition wu prior to the descent of the cold oceanic litho- 
sphere. 

One of tile most prominent features in plate convergent ereas is the island 
arc volcanism. Unfortunately, due to the coarse grid system of the present 
models, a detailed study of island arc magma generation is not possible. 
Since both surface topography and gravity in these areas are dominated by 
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20 0 MT 



Fig. 1 1 . Santa at Fig. 7 axcapt alab aubduction ia initiated at point A in Fig. 1 . 









iilind arcs, the failure to create island arc volcanoes prohibits a meaningful 
{Mresentation of the surface topography and gravity anomalies predicted by 
the present models. However, the present model is able to confirm that the 
induced flow in the wedge area is indeed a strong convection flow, as dis- 
cussed earlier. Even though by itself alone, the induced flow may not be 
sufficient to initiate back arc spreading (Toksdz and Hsui, 1978), it remains 
a plausible mechanism nevertheless. 

SUM.MARY 

A thermal regime of a slab subducting into the mmtle and its surrounding 
regioiis has been calculated taking into account the effect of convection with 
temperature, pressure dependent viscosity. Results indicate that the subdue- 
tion process is primarily a cooling process. The cold descending lithosphere 
and the induced flow are the dominating factors in determining the thermal 
structure of the region at any given time. The temperature within the slab 
and in its immediate environment is insensitive to the mantle condition prior 
to slab subduction. In order to maintain a high temperature region beneath 
the inter-arc basins, as suggested by the high heat flow and the high attenua- 
tion of seismic waves, our model indicates that viscous sheiur heating within 
the mantle is probably an important factor. Without this, the wedge area is 
likely to be cooled rapidly by the descending cold lithosphere. Our model 
of the thermal structure of a downgoing slab is in good agreement with other 
slab models and it represents an average of these calculations. Finally, the 
induced flow within the wedge area is indeed a prominent feature of the sub- 
duction process. Its contribution to the mechanism of inter-arc spreading 
remains viable. 

The primary purpose of the present model was to understand the thermal 
evolution of a slab and its surrounding mantle over a broad region. Because 
of the coarse grid system used for computational efficiency, detailed features 
such as magma generation and island arc volcanism cannot be studied. In 
order to have models which can explain detailed surface features, further 
work with improved computational techniques that yield finer resolutions 
must be carried out. In the mean time, the models presented here provide a 
better general understanding of the subduction process and th^y represent an 
advancement towards a fully dynamic and more realistic model of the sub- 
duction regions. 
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Subduction of continental lithosphere 
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ABSTRACT 

Two afftcts that contributa to th* aubduction of continental Utho* 
sphere arc the negative buoyancy of the relatively cold mantle part of 
continental lithosphere and the pull of a downgoing slab of oceanic litho* 
sphere on continental lithosphere trailing behind it. The mantle part of 
the lithosphere could continuously subduct about 10 km of continental 
crust, if the upper and lower crust could be detached from one another. 
This estimate, however, could be in error by as much as a factor of three, 
given the range of plausible values for the requisite parameters. The second 
effect is more difficult to estimate because of our ignorance of what pro* 
portion of the gravitational force acting on the downgoing slab is trans- 
mitted to the surface lithosphere. For various assumptions, subducted 
oceanic lithosphere could be expected to pull a short (tens of kilometres) 
or long (hundreds of kilometres) length of intact continental crust into 
the asthenosphere. Although we favor subduction of only a short length 
of intact continental crust, we cannot prove that a large amount of crust 
is not subducted. These calculations are for continental margins that meet 
subduction zones flush. If continental crust can remain intact, peninsulas 
and n^crocontinents might be subducted completely. 


INTRODUCTION 

Because of its low density, continental crust has traditionally 
been assumed to be too light to be subducted. Even with the recog- 
nition of the important role of the lithosphere in subduction, 
eanh scientists have assumed that the large thickness of continental 
crust provides too much buoyancy to continental lithosphere for 
it to be subducted (for example, backs and others, 1968; McKenzie, 
1969). Yet geologic evidence indicates large-scale crustal shorten- 
ing by thrust faulting in most mountain belu, particularly where 
continents collide. Consequently, some subduction of continental 
crust must take place, and in general, geologic evidence pbcesonly 
lower bounds on the amount of thrust displacement and therefore 
of subduction of lower crust. For instance, Cansser (1966) esti- 
mated that the total displacement on the major thrust faults in the 
Himalayas is at least 300 km, excluding displacement at the Indus 
suture '.one, where thousands of kilometres of oceanic lithosphere 
probably were subducted. Powell and Conaghan (197.')), however, 
proposed as much as 1 ,000 km of slip on one of the Himalayan 
faults, the Main Central thrust. In the Norwegian Caledonides, 
Gee (1975) inferred at least SOO km of shortening distributed on 
several faults. These estimates do not necessarily require a large 
amount of subduction of the entire intact continental crust to 
depths of hundreds of kilometres in the asthenosphere. In all 



cases the thrust faulu are exposed only within the crust and could 
be interpreted as detachment surfaces between the upper and lower 
crust. Nevertheless, a substantial amount of crustal material must 
be subducted if only to depths of 70 km, to account for the 
observed shortening. 

A basic question is. If the buoyancy of continental crust 
does inhibit subduction of continental lithosphere, then how much 
crust can be subducted before convergence is halted? Resolution 
of this question will probably require further detailed geologic 
mapping of thrust belu and may require new techniques. At 
present, however, we can examine limits on how much crust could 
be subducted given basic assumptions about the physical proc- 
esses and requisite parameters involved. 

Two phenomena contribute to the subduction of continental 
lithosphere. First, the mantle part of the lithosphere is cold and 
therefore dense enough that its negative buoyancy would pull the 
bottom part of the crust down if the top of the crust could be 
detached from it. Second, the gravitational body force acting on 
the subducted oceanic lithosphere may also exert force on the 
leading edge of a continent following the oceanic lithosphere 
into the subduction zone. 


geology. 7. 0 . 58-62. JANUARY 1979 



- A* 






GRAVITATIONAL STAIILITY 
OP CONTINENTAL UTHOSPHIRl 

Oceanic Uthoiphcri plungct into the aMhtmoeplMrt because 
it is gravitationaliy unstable. It it colder and tnarafore more dense 
than the atthenospherc would be if brought adiidMikally to the 
tame depth (pressure). Although the crust it lighter than the mantle 
and buoys the lithosphere up, evidently oceanic crust it too thin 
to prevent tubduction. The effects of th«e two buoyancy forces 
can be expressed simply u the differmce bmweea two terms, one 
for the buoyancy of the crust and the other for the negative 
buoyancy of the mantle part of the lithosphere (Fig. 1; McK«i^, 
196^. For the crust, the buoyancy force per unit area inhibiting 
subduoion is (,p„ - Pe)|h, where and are the densities of 
the mantle and crust, g is gravitation^ acMteation, and A is the 
thickness of the crust. For the mantle part of the Uihosphme, the 
difference in the dnuity between material with tamperuure 7{z), 
and that with temperature Tgit), following the adiabmic gradlmt, 
is simply (Tte) - Tg{t)\, where a is the coefficient of 
thermal expansion and < is the depth. 

The temperature in the uthenotphere presumably dosely 
follows the adiabatic gradient (BO.S'/km). Beatuse of the ti^l 
amount of radioactivity in the mantle, the geotherm in the litho- 
sphere is essentially linear (Fig. 1). If the base of the lithosphere 
is defined by a particular isotherm, Tg, then in oceanic lithosphm, 

nt)*^x, ( 1 ) 

where a is the thickneu of the lithosphere. As the icmperaturit in 
only the mantle part of the Uthosphere is important, a mors pre- 
cise description of the temperature profile across it is 

nz)-{Tg~Tc)^*Te, ( 2 ) 


depth of the Moho (Fig. t). For conUnental lithosphere, equation 
2 to tUghUy batter t^ equation I. because radioactivity in the 
continanul am makes an important contribution to the geo- 
dsma. The spproxlmatton of the gaothmra in the mantle as two 
swemislly straight lines— one due to conduoimi through the 
itohosphen and the other u an adiabat maintaiiMd by coo- 
vaettom—to, of courts, a simpUfleation. More rMitotic gaotherms, 
dtoeussed by Parsons and McKenzio (I97g), however, are not 
srarrantad here. 

The negative buoyancy force per unit area due to the low 
temperature of the lithosphere to Pmagf'h ~ HOIde. 
This can be written directly by reealung the formula for a triangle 
(see Pig. I): -HAlKe - A). Because the adia- 

batic gradient it to small, Tg(g) to nearly constant, so that 7^(A) 
■ V«) ■ tmnperature of the asthenotphere. 

For an insubility to exist, the mantle part of the lithosphere 
must be haaviw than the crust, or 

- Pc)F* - I Tgia) - r(A)l (e - A)< 0. (3) 


Reasonable values for the various parameters are 
p„ ~ Pg m 3.33 - 2.1 ■ 0.39 * 0.1 g/cm* 

P„ m 3.39 w 0.09 g/cm’ 
a - 3 * I K I0-* ‘C** 

Tg m 1200 * 190 *C 
7Xhg) ■ 900 * 100 *C (for continents) 
e - Ac * >S0 w 90 km (for continenu) 
e - Aq ■ 100 km (for old ocean fltoor) 

For oceanic lithosphere, where /^ « a ■ 100 km, ■ 0 
*C, and pg m 2.9 g/em*, and insubility exitu when 


, Pm® Ifa(a)- HAo)) (a - hg) 
® 2(Pm - Po) 


< 13.4 km. 


where Tg is the temperature at the base of the crust and A is the 

iitmiitruiu 
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Figure I. Geothoims and gravitational body (orcs. In Uihosplicre, 
heat is presumably conducted, and icmpcraiure gradient is essentially 
I lineal. Beneath liUiospliere in aslhenosphcie, heat isconvected, and 
\padicni is adiabatic. Because lithosphere is colder than asthenosphere, 
it is more dense, and a body force due to gravity acts on ib This negative 
buoyancy force is opposed by buoyancy of crust, which is le» dense than 
maniic. Integral over depth of curves on righb for oceanic and conti- 
nental Uthosphere. gives force per unit area acting on lithosphere. 


This condition to nsily fulfilled in oceanic regions. 

Expression 3 allows an estimate of what fraction of the lower 
continental crust, hp would be unstable if the upper part could be 
detached from it. Using the values given above, hg < lO km, ap- 
proximately 309ii of the crust. This estimate is very uncertain, 
however. Note that if we choose the extreme values given above, 
Af < 2 km or A, < 29 km. In the former ease, a small fraction 
of the crust would be carried down, whereas in the latter nearly 
the entire continental lithosphere would be unstable! These ex- 
treme estimates arise from a conspiracy of extreme values for 
each of the parameters given above, and therefore both estimates 
are unlikely to apply in general. The calculations do show, how- 
ever, that we cannot eliminate the possibility of subducting a large 
fraction of continental crust, if it could be detached from the 
upper part. 

GRAVITATIONAL FORCE ACTINO ON DOWNCOINC 
SLAB OF OCEANIC LITHOSPHERE 

If a continental mass is part of an oceanic plate being 
consumed at a subduction zone, gravitational body forces acting 
on the excess mass in the slab might be able to pull continental 
crust into the asthenosphere (see Fig. 2). Given plausible param- 
eters that govern the various forces, we can place an upper bound 
on what length (d) of continenul crust could be subducted. There 
are basically two questions here: What is the gravitational force 
acting on the slab? How much of it actually pulls on the part 
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Figun 2. Gravittdonal fom Kting on ooMnfe UthM^twn md 
pulUng eoniinmtil litfioiphtrt into utfitnotphm Fom it oppo«d 
by btioytiiey of eontinmtil emit of , 


of th« lithoipherc at cht surface? Given our inaWlity to answer 
the latter question, the uncertainty in the estimate of the amount 
of continental crust that can be subducted is very large. Never- 
theless, we present the resuiu for an assortment of usumptions 
and then state our prejudices about them. 

The gravitational body force (pw unit length of subduetion 
gone) acting on the slab is due to the temperature of the slab being 
lower than that of the surrounding asthenosphere, and it an be 
estimated from calculations of the temperature in the slab, which 
depends primarily upon the subduetion rate and the plate thickness 
(for example, McKenzie, 1969). Numerical caiculations'show that 
radioactivity in the oceanic crust and mantle does not contribute 
much to the temperature in the slab (Minear and ToksOz, 1970). 
The role of frictional heating remains controversial, but since we 
seek an upper bound on the amount of crust subducted, here we 
ignore its contributions both as a heat source and as a retarding 
force. Once a convergence stopsr friction ceases to be either a 
heat source or a retarding force. We discuss phase changm beirw. 
The principal source that heats the slab is conduction from the 
surrounding asthenosphere (McKenzie, 1969; Minear and 
ToksOz, 1970). 

McKenzie (1969) derived an analytic solution for the 
temperature in the slab and calcuiated the gravit.ttional 'orce 
acting on the slab due to iu excess mass. This cal.uiation is given 
in the Appendix 1, but it is extended to include different flnite 
depths of penetration of the slab, and it is evaluated using more 
reasonable parameters for the plate thickness and surrounding 
temperatures than McKenzie used. From equation A5 and the 
parameters in the appendix, we obtain forces (per unit length 
along the arc) of 1.8 x 10'* dyne/cm for a slab extending to a 
depth of 350 km (than is, 500 km long and dipping at 45*) and of 
3.2 X 10" dyne/k-m for a slab extending to a depth of 7cio km 
(that is, 1 ,000 km long and dipping at 45*). For comparison, we 
note that Fitch (1977) integrated the force per unit length calcu- 
lated by Schubert and others (1975) and obtained a maximum 
streu of 3.2 kb acting on a plate 100 km thick extending to 700 
km depth, corresponding to a force per unit length of 3.2 x 10'* 
dyne/cm. A crude estimate of the excess mass calculated by 
ToksOz and othen (1973) yields a similar value. 

The principal sources of error in ail of these calculations are 
the values of the thermal expansion, a, and the thicknesses of the 

3 ^ 2 -/, 


ptaia, (boeiust tha force p» unit length is proportional to Oq*). 
lemise of the uncertainties in the parameten inserted into these 
formulas, the resulting calculatioiu could be in error by u 
much as a factor of two. 

Another possible force is the excess mass in the slab due 
to the possible occunence of the oUviiM-spiael phase change at 
a shallower depth in the slab than in the surrounding anhsoo- 
rphere (for example, Richter, 1973: Sehubm and othm, 1975). 
Rich (1977) integrated Schubm and others* calculation to obtain 
an additional force of 1.9 x 10" dyne/km. There it a Mntrary 
suggestion, however, that the reaction h inhibited by iu kinetics 
In the cold central part of the slab the reaction would be to slow 
that this pan of the slab would pmetrats well below the equilib- 
rium depth before the reaction would o«ur (Sung and Bums, 
1976a, 1976b). In this cate, it wculd not add a force pulling the 
slab down; instead, it would oppose subduetion (Sung and 
Bums, 1976a, 1976b). 

It it clear the the uncertainty in tbs gravitational body 
force acting on a downgoing slab it laige. la addition, ih«e are 
few dau contuaining how much of this force it transited to 
the lurflcial pan of the lithosphere. Some smms man in the slab 
will be tupponed by stress in the tuitounding atthenotphm. 

For instance, the fault-plane toiutiont of deep earthquakes sug- 
gest that the exceu mass of downgoing tlabt below about 330 km 
does not pull the shallow pan down (Isacks and Molnar, 1969, 
1971). Perhaps the force acting only on the pan shallower than 
350 km pulls the surflcial lithosphere. In fact, the pouibUity exitu 
that the downgoing slabs do not exen a substantial force on the 
surface lithosphere. Kanamori (1971) suggested that at least mme 
slabs break off at the trench and sink freely into the astheno- 
sphere. As the rate that the lithsphere sinks into the anheno- 
sphere is comparable to the rate at which it would sink into a 
fluid with the vUcosity of the astheniuphere, Richter (1977) 
inferred a small force exened by the slab on the rest of the plate. 
Richter and McKenzie (1978) suggestnl that the downgoing slab 
exeru a maximum stress of only 100 bar on the remaining surflcial 
lithosphere. For a plate 100 km thick, this leads to a force per 
unit length of 10" dyne/cm. Like the others, this estimate is 
surely uncertain by a factor of two of more. 

Thus, the forre per unit length can be very different 
depending upon the assumptions made and can be summarized 
as follows: 

F w l.C X 10'* dyne/cm: little direct force applied by the 
sinking slab to the surface plate 

F a 1.8 X 10'* dyne/cm: pull from a slab to a depth of 
350 km 

F a 3.2 X 10" dyne/em: puil from a slab to a depth of 
700 km 

F a 5.1 X 10'* dyne/cm: pull from a slab to a depth of 
700 km with excess mass in an elevated olivine-spinel 
phase change 

The variation in the values underscores our uncertainty in 
the proceues occurring. 

LENGTH OF CONTINENTAL CRUST SUBDUCTED 

The pull of the slab is counteracted by the buoyancy of 
the light continental crust. This force per unit length is (p„ 
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- (• iht lM|ih of (te co nt laww i i cnitt Uuu 

b subducttd (Fig. 2). Tbt mintU pvt of Um c o ii Bcm il Utho* 
iphcrt will coatribiKt to tht nibductkm by pulliag down • frae* 
linn of Uw oiitt, hg. Tbt hatvy ocouic slab b tbwtfoft rwiuad 
by (p„ • Pe)t(t>e ** V^’ eoMiofaiul Uibofpbtra 

ibould CMM whtn tht rorct pulUng it dowii b biiaiicad by buoy* 
incy of tht cootincnul enut. Wt sack ■ vthit of 

ii m ^ , 

(fim “ Pc)f(*c • 

tbt Itngth of tht comiatnt nibduotd. Uiiiif (p„ -> p^) ■ 0.5S 
g/cm*. g ■ to* cm/I*, tad (b( ^ ■ 3S • 10 ■ 25 km, wt 

obtain for ^ 1.0 x 10** dyM/em: d ■ 7 km; ^ ■ 1.1 x 10'* 

dynt/em: d ■ 131 km; F ■ 3.2 x 10“ dynt/cm: d •Til km; 
and F ■ S.l x 10'* dynt/cm: d ■ 2M km. Noit that an error 
of 0.1 g/cm* in (p„ • pf) contributM a t9dk error in tht Iragth 
of the enut subducted. Recall abo that the uDctnaimy in hg 
calculated in tht previous section w ,.\ld nearly allow ttaady itate 
lubduction of continental liihoephere. If ■ 20 km inst^ of 
10 km, tht amount of mut to ^ subducted would be about 
40^ greater than that calculated above. As in the previous cal- 
culation, the 33ak uncertainty of a in the paramettn listed above 
leads to another 33F« uncertainty in d. Each of these estimates 
for the amount of continental crust subducted b therefore uncer- 
tain by at least a factor of two, but notice that the uncertainty 
due to our ignorance of tht processes involved leads to a much 
larger uncertainty in d. 

These calculations arc ail based on the implicit assumption 
that the continent meets the subduction zone flush and thu 
throughout the length of the subduction zone there is continental 
material. Such a situation b highly unlikely, for surely penituulas 
will meet the subduction zone first. Moreover, continental frag- 
menu need not have the same width as the length of the subdue* 
tion zone. The buoyancy of a long peninsula or of a small con* 
dnental fragment, such u the Ssycheilerlslands, b surely negligible 
compared with the negaUvt buoyanee of the oceanic lithosphere, 
subduction of them ought to be coniplete. In the case of a conti- 
nent half as long as the subduction zone, the amount of continent 
•iiat could be subducted should be twice the values given above. 

We think that except for such casm, howevei , only very 
short lengths (d - SO km) of ivaef continental cruet can be sub- 
ducted. The gravity anomalies across bland arcs (Kogan. 1973; 
Watts and Talwani, 1973) do not show evidence for the excess 
mass of the elevated olivine-spinel phase change. Because the 
fault-plane solutions of earthquakes below 330 km imply a 
resistance to the sinking of the slab, this lower ('art probably 
does not exert a puil on upper parts. Therefore, the shortest 
length (131 km) seems the most plausible of the three largest 
values given above. 

Beyond this, the evidence that large events beneath the 
trenches seem to rupture through the entire lithosphere (Abe, 
1972; Kanamori, 1971; Stewart, 1978) and the theoretical argu- 
ments of Richter (1977) and Richter and McKenzie (i978) suggest 
that only a small fraction of the negative buoyancy of the sinking 
slab is transferred directly to the surface. Moreover, the outer 
(topographic rises seaward of deep-sea trenches do not seem to 
Require large forces or moments that would be applied by the 
weight of the downgoing slab. Using laboratory data to constrain 
the strength of the lithosphere, C. Goeue and B. Evans (in prep.) 
have concluded that a bending moment per unit length of about 


10** dyiM could eauM a flaxurt of tht Uthosphere and ihareby ac- 
count for the obterved topography. If ihk momont wweeauMd by 
tht ne^oivt bttoyiney force par unit langih of tbt downgoing slab. 
tppUad 100 km landward of tht trmeb tab. a force ptr unit 
Iti^ of 10** dynt/em would bt adaquatt. A much larger force 
could deflect the oceanic litboepime much more than b observed 
If thb logic b correct, thm only a short length of intact con- 
tinmital crust would b« subducted. 

Subd.uction of a short length of intact crust b not necessarily 
incomimeibie with the large oiittal shortaoing obswved in oro* 
genic bdti. The geoiogie evidence impUai detachment of the 
upper and lower crust along shaliow-d^plng thrust fauhs. It the 
upper crust can be detached, (he negative buoyancy of the mantle 
part of continental lithosphere b likely to be large enough to 
sustain subduction of it and the lower crust. Moreover, if the 
crust has a low melting tmnperatura, it might melt before pene- 
trating deep into the asthenosphere and thereby fkirther redwe 
the resbtanca to subduction. 


A DIGRESSION ON PRECAMBRIAN BLUESCHI8TS 

Biueschbu arc most common in Mesozoic and Canozoic 
terrenes; there seem to be no Prerambrian examples (de Roever, 
1936; Ernst, 1972). Greater heat gmieration and higher ^thermal 
gradienu in Prccambrian time are thought to have prevented the 
high-pressure-low-tcmperature tubility field for biueediisu at 
that time (de Roever, 1996; Ernst. 1911). Anothw possibility b 
that ftecambrian lithosphere was too thin and too buoyant to bo 
able to pull cold crustal matmiai down to the depths required 
for formation of biueschbu. 


SUMMARY AND CONCLUSIONS 

Two factors contribute to subduction of continenul 
lithosphere: the negative buoyancy of the mantle part of the 
lithosphere and the force exerted on the surface lithosphere by 
a downgoing slab of oceanic lithosphere. If there were a mechan- 
ism for detaching the upper crust from the lithosphere, the nega- 
tive buoyancy of the mantle part would make it and the lower 
ciust graviutionally unstable. The mantle lithosphere could 
subduct about 10 km of crust with it, but thb estimate of 10 km 
is very uncertain and could be nearly three times larger if the 
values of the requisite parameters conspbed. Gravity acting on 
the excess mass of a downgoing slab of lithosphere exerts a 
downward force on the surficial lithosphere. This force should 
be able to puil some continental lithosphere into the astheno- 
sphere. If the crust remains intact, peninsulas and microcontin- 
enu might be subducted completely. The greatest uncertainty in 
the calculation of this force arises not in estimating the excess 
mass but in deciding how much of its gravitational body force 
b transmitted to the surface. If all of the force were transmitted, 
then hundreds of kilometres of continent could be subducted, as 
noted also by Bird and others (1973). If only a small fraction 
were transmitted, then only a negligible length of intact contin- 
ental crust could be subducted. Although we think that the latter 
b the case, at present we cannot demonstrate it conclusively. If 
a substantial amount of continental crust can be subducted, (his 
clearly will cause pronounced chemical heterogeneity of the 
mantle and will profoundly affect the crustal evolution of 
subduction zones such as in Tibet. 
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Thil dtiivation waiyivcn by McKtnett (1969), but he let •**, M that tht 
turn rn equation AS vanithev The number, in, uf termi uied in turn dependi 
upon the depth ut penetraiiun ui the Hab. but lor fewer than me 

lermt are adequate, and tor x,, lb 3t. only the first term is important. 
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SOHB CONSTRAINTS ON THE SHEAR NAVE VELOCITY 
STRUCTURE UNDER THE TIBETAN PLATEAU 

W.P. Omn 

P. l4oinar (both att Dopartiiwnt of Earth and Planatary Sclancaa 

Haaaachuaatta Znatituta of Tachnolcgy* CambrldgOf Maaa. 02139. 

Wa hava eomblnad tha obaarvatlons of aurfaca wava dla- 
paralon acroaa tha Tlbatan plataau and tha dlffarantial traval 
tlma raalduala batwaan P and S wavaa for aarthquakaa in Tibat 
to placa bounda on tha poaaibla ahaar wava valoclty atructura 
banaath Tibat. Tha P wava valocity atructura ia ralativaly 
unconatrainad and waa aaaunad to ba tha aama aa Dziawonaki'a 
valocity modal S2 for tha ahiald ragiona axcapt for tha in- 
cluaion of a poaaibly thick cruat. Surfaca wava diaparaion 
data auggaat a ranga of poaaibla cruatal thicknaaa from 55 km 
to 70 km or avan 85 km with upper mantle ahaar wave valocitiaa 
of 4.4 km/aec, 4.7 km/aacr and 4.8 km/aac raapactivaly. Tha 
travel time reaiduala indicate the ahear wave velocity beneath 
Tibet ia, on the average, 5 to 6% leaa than that beneath ahielda, 
asauming that all the differencea in velocity atructura between 
Tibet and ahielda were confined to the upper 240 km. If the 
ahield model S2 ia adjuated to have the same cruatal thickneaa 
aa the Tibetan modela, the Tibetan modela will have average 
ahear wave velocitiea 4 to 5% leaa than that of S2 in the 
depth range of 100 km to 240 km. Thia correaponda to an 
average ahear wave velocity of about 4.3 km/aec in thia depth 
range, and might be an indication of elevated termperature. If 


• strong low-Q structuro •xists undor Tlbot* thon bocaus* 
of velocity dltperalon, the average velocities should be 
even lower. 


(Presented at the Pall Annual Meetings of the American Geo- 
physical Union, December, 1978.) 


SUBDUCTZON INDUCED CONVECTION AND 


ISZAND-ARC CONVECTION 


Albyt T. HBUi 

>4. t^afi tokiOz (both att Dapartmant of Earth and Planatary Sciancaa^ 

Maaaachuaatta Inatltuta of Tachnology, Cambridga* MA 02139 
Bruca D. Marsh (Dpaartnant of Earth and Planatary Sciancas, Johns 

Hopkins Unlvarslty, Baltlmora, MD 21218) 

Subductlon of lithospharic slabs inducas co.ivactlon 
in tha ovar lying asthanosphara. This convaction strongly 
influancas tha tharmal ragima of tha subductlon ragion, and 
has implications to tha formation of back-arc basins and 
island-arc volcanism. In this study, wa carry out a datailad 
calculation on tha thermal structure of tha wedge area between 
tha subducting slab and the overlying lithosphere, utilizing 
a finite difference numerical technique. Results indicate 
that the area above the descending lithosphere will generally 
be of higher temperatures if a low viscosity asthanosphare 
exists. This is mainly because the low viscosity astheno- 
sphere will deform in response to the shear stresses applied 
at the boundary by the descending lithosphere. As a result, 
and induced convective current is established above the slab. 

This current will steadily feed the wedge corner with warmer 
material from tha deeper interior. The 'cooling effect* of 
the cold slab can only penetrate to a certain skin depth which 
is related to the subductlon speed. One the other hand, at 
the absence of convection, the wedge area will continuously 
be cooled by the descending lithosphere. Therefore, from 


th« thsmal point of vlow alono* it is apparsnt that island 
are volcanism is more likely to take place if there exists# 
above the subducting lithosphere# a low viscosity region where 
convective currents can be induced. This is consistent 
with the conclusions derived from seismic observations. 


(Presented at the Fall Annual Meetings of the American Geophy~ 
sical Union# Oeconber# 1978.) 


